arXiv: 1506.0517lv3 [q-bio.QM] 21 Jan 2016 


On the Origins and Control of Community 
Types in the Human Microbiome 


Travis E. Gibson^’*, Amir Bashan^’*, Hong-Tai Cao^, 
Scott T. Weiss^, and Yang-Yu 


1 Channing Division of Network Medicine, Brigham and Women’s Hospital, 
Harvard Medical School, Boston, Massachusetts 02115, USA 

2 Department of Electrical Engineering, University of Southern California, 
Los Angeles, California 90089, USA 

3 Center for Cancer Systems Biology, Dana-Earber Cancer Institute, 

Boston, Massachusetts 02115, USA 

* These authors contributed equally 

f Corresponding Author: yyl@channing.harvard.edu 


Microbiome-based stratification of healthy individuals into compositional 
categories, referred to as “community types”, holds promise for drastically 
improving personalized medicine. Despite this potential, the existence of com¬ 
munity types and the degree of their distinctness have been highly debated. 
Here we adopted a dynamic systems approach and found that heterogeneity 
in the interspecific interactions or the presence of strongly interacting species 
is sufficient to explain community types, independent of the topology of the 
underlying ecological network. By controlling the presence or absence of these 
strongly interacting species we can steer the microbial ecosystem to any de¬ 
sired community type. This open-loop control strategy still holds even when 
the community types are not distinct but appear as dense regions within a 
continuous gradient. This finding can be used to develop viable therapeutic 
strategies for shifting the microbial composition to a healthy configuration. 



2 



Contents 


1 Main Text 5 

1 Introduction. 5 

2 Dynamic Model. 6 

3 Metacommunity and Local Communities. 7 

4 Origins of Community Types. 8 

5 Discussion. 10 

2 Figures 13 

3 Materials and Methods 19 

4 Supplementary Figures 23 

5 Supplementary Text 39 

1 Introduction. 39 

2 Notation. 39 

3 Random Variables and Random Matrices . 40 

3.1 Primer on Random Variables. 40 

3.2 Random Network Models. 42 

3.3 Spectrum of Random Matrices. 43 

4 Stability. 44 

4.1 Primer on Stability. 44 

4.2 Stability of Generalized Lotka Volterra Dynamics. 46 

4.3 Stability in the presence of new species . 47 

4.4 Robustness to disturbances. 49 

4.5 Diagonal Stability of Random Matrices . 58 

5 Clustering Analysis and Ordination Methods. 59 

5.1 Distance and Metrics. 59 

5.2 /c-Medoids. 59 

5.3 Silhouette Value. 59 

5.4 Variance Ratio Criterion and the Calihski-Harabasz Index . 60 

5.5 Principle Coordinate Analysis. 60 

5.6 Principle Component Analysis . 61 

6 Modeling . 61 

6.1 Universal Model. 61 

6.2 Mutiple Set Model. 63 

7 Simulation Results and Analysis . 63 

7.1 Multiple Set Models. 63 

7.2 Universal Model. 63 

Bibliography 97 


3 
































4 


CONTENTS 



Chapter 1 


Main Text 


1 Introduction 

Rather than simple passengers in and on our bodies, commensal microorganisms have 
been shown to play key roles in our physiology and in the evolution of several chronic 
diseases [25,81]. Many scientific advances have been made through the work of large- 
scale, consortium-driven metagenomic projects, such as the Metagenomics of the Human 
Intestinal Tract (MetaHIT) [71] and the Human Microbiome Project (HMP) [94,95]. In 
particular, the HMP has analyzed the largest cohort and set of distinct, clinically relevant 
body habitats to date, in order to characterize the ecology of human-associated microbial 
communities [95]. These results thus delineate the range of structural and functional 
configurations normal in the microbial communities of a healthy population, enabling 
future characterization of the translational applications of the human microbiome. 

A recent study proposed that a healthy gut microbiome falls within one of three distinct 
community types, which the authors coined as “enterotypes” [3]. More specifically, the 
authors calculated the relative abundance profiles of microbiota at the genus level and then 
performed standard cluster analysis, finding three distinct clusters (enterotypes). Each 
enterotype is a dominated by a particular genus {Bacteroides, Prevotella, or Ruminococcus) 
but not affected by gender, age, body mass index, or nationality of the host. These 
results suggest that enterotyping could be an efficient way to stratify healthy human 
individuals. The development of personalized microbiome-based therapies would then 
simplify to shifting an unhealthy microbiome to one of the distinct healthy configurations. 

A met a-analysis, however, suggested that enterotypes, or in general community types, 
could be an artifact of the small sample size in [3] and what one should expect is a 
continuous gradient with dense regions rather than distinct clusters [52]. The level of 
discreteness or continuity of the community types remains unclear. Interestingly, samples 
in the dense regions of this gradient are either highly abundant or deficient in Bacteroides 
[52], indicating that community types could still emerge as the dense regions within a 
continuous gradient. Indeed, some recent work actually supports the notion of distinct 
community types [22,41,79,102,107]. 

We still lack consensus on the nature and origins of community types [4,14,32,45,51]. 
In principle the presence of community types can be explained by several different mech¬ 
anisms. First, there may be true multi-stability, i.e. multiple stable states with all mi¬ 
crobial species present in the same environment [57]. Those stable states are interior 
equilibrium points of the corresponding ecological dynamic system. Although this type 
of multi-stability has been well discussed in macro-ecological systems [17], its detection in 
host-associated microbial communities is rather difficult and has not been demonstrated 
experimentally, partially because any subtle differences in the environment can drive those 
microbial communities [32]. Second, there may be strong host heterogeneity, leading to 
host-specific microbial dynamics (parameterized by host-specific intra- and inter-species 
interactions). If those interactions, which serve as parameters of the host-associated micro¬ 
bial ecosystems, can be classified into distinct groups, then we can numerically demonstrate 
that distinct community types will naturally emerge (Supplementary Text Sec. 6.2 and 
7.1). Yet, the presence of classifiable microbial dynamics has not been experimentally 
detected, presumably due to the lack of high-quality time-series data for a large number of 
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subjects. Moreover, the overwhelming success of Fecal Microhiota Transplantation (FMT) 
in treating recurrent Clostridium Difficile Infection (rCDI) actually implies that difference 
in conditions between individuals are unlikely the cause of community types [98,105,106]. 

Here we proposed a simple mechanism, without assuming multi-stability or host hetero¬ 
geneity, to explain the origins of community types. In particular, using a dynamic systems 
approach, we studied compositional shift as a function of species collection and demon¬ 
strated that with heterogeneous interspecific interactions, a phenomenon often observed 
in macroecology [26,27,73], community types can naturally emerge. Interestingly, this 
result is independent of the topology of the underlying ecological network. To our knowl¬ 
edge, this is the first quantitative attempt to explore the analytical basis of community 
types. Furthermore, community types, even when they weakly exist, can be manipu¬ 
lated efficiently by controlling the Strongly Interacting Species (SISs) only.^ This provides 
theoretical justification for translational applications of the human microbiome. 

2 Dynamic Model 

The human microbiome is a complex and dynamic ecosystem [35]. When modeling a 
dynamic system we should first decide how complex the model needs to be so as to capture 
the phenomenon of interest. A detailed model of the intestinal microbiome would include 
mechanistic interactions among cells, spatial structure of the human intestinal tract, as 
well as host-microbiome interactions [10,68,83,97]. That level of detail however is not 
necessary for this study, because we are primarily interested in exploring the impact that 
any given species has on the abundance of other species. To achieve that, a population 
dynamics model such as the canonical Generalized Lotka- Volterra (GLV) model is sufficient 
[32,76]. Indeed, GLV dynamics leveraging current metagenome data has already been used 
for predictive modeling of the intestinal microbiota [33,62,86]. 

Gonsider a collection of n species in a habitat with the population of species i at time 
t denoted as Xi{t). The GLV model assumes that the species populations follow a set of 
ordinary differential equations 

n 

Xi(t) = riXi{t) + Xi{t)'^aij Xj(t), z = 1,... ,n (1) 

i=i 

where () = ^ 0 - Here is the growth rate of species z, aij (when i ^ j) accounts 
for the impact that species j has on the population change of species z, and the terms 
aiix‘i are adopted from Verhulst’s logistic growth model [37]. By collecting the individual 
populations Xi(t) into a state vector x{t) — [xi(t), ••• ^Xn(t)Y ^ Equation (1) can be 
represented in the compact form 

x{t) = diag(x(t)) (r + Ax(t)) , (2) 

where r = [ri, • • • , XnY is a column vector of the growth rates, A = ffiij) is the interspecific 
interaction matrix, and diag generates a diagonal matrix from a vector. Hereafter we drop 
the explicit time dependence of x. 

Next we discuss the notion of fixed point, or equivalently steady state, in the GLV 
dynamics. This notion is important in the context of the human microbiome, as the mea¬ 
surements taken of the relative abundance of intestinal microbiota in the aforementioned 
studies typically represent steady behavior [3,95]. In other words, the intestinal microbiota 
is a relatively resilient ecosystem [61,80], and until the next large perturbation (e.g. antibi¬ 
otic administration or dramatic change in diet) is introduced, the system remains stable 
for months and possibly even years [13,20,31]. The fixed points of system (2) are those 
solutions X that satisfy x = 0. The solution x = 0 (i.e. all species have zero abundance) is 
a trivial steady state. The set of non-trivial steady states contains those solutions x* such 
that r + Ax* = 0. When the matrix A is invertible, it follows that the non-trivial steady 
state X* = —A~^r is unique [38]. 

Our study ultimately investigated the impact that different collections of microbial 
species have on their steady state abundances. In Box I we presented a detailed analysis 

^In this paper we use the term species in the general context of ecology, i.e. a set of organisms 
adapted to a particular set of resources in the environment, rather than the lowest taxonomic rank. One 
could think of organizing microbes by genus or operational taxonomical units as well. 
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showing that if we introduce a new species into the ecosystem in (2), the shift of the steady 
state is proportional to the interaction strengths between the newly introduced species and 
the previously existing ones. Similarly, if two communities with the same dynamics differ 
by only one species, then it is the interaction strength of that species with regard to the rest 
of the community that dictates how far apart the steady states of the two communities will 
be. This analytical result indicates that heterogeneity of interspecific interactions could 
lead to the clustering of steady states, and hence the emergence of community types. 

To systematically investigate how changes in species collection affect the steady state 
shift in the GLV dynamics, we assumed that two microbial species will interact in the 
same fashion regardless of the host. Otherwise, if the interactions are host specific and the 
dynamics are classifiable, we can show that distinct community types will emerge almost 
trivially (Supplementary Text Sec. 6.2 and 7.1). 

3 Metacommunity and Local Communities 

Consider a universal species pool, also referred to as a metacommunity [18], indexed by 
a set of integers S = {1, ..., n}, annxn matrix A representing all possible pairwise 
interactions between species, and a vector r of size n containing the growth rates for all 
the n species. The global parameters for the metacommunity are completely defined by 
the triple (S,A,r). We consider q Local Communities (LCs), defined by sets that 
are subsets of S, denoting the species present in LCu with u = 1,... ,q. This modeling 
procedure is inspired by the fact that alternative community assembly scenarios could 
give rise to the compositional variations observed in the human microbiome [18]. These 
LCs represent microbial communities in the same body site across different subjects. For 
simplicity, we assume that each LC contains only p species (p < n), randomly selected 
from the metacommunity. 

The CLV dynamics for each LC is given by 

LC^y : = diag T A ^'^^, (3) 

where the LC specific interaction matrix and growth vector are defined as 
and respectively. That is, is obtained from A by only taking the rows 

and columns of A that are contained in the set A similar procedure is performed in 
order to obtain Finally for each x^^^ there is a corresponding G MU that has the 
abundances for species of LCu in the context of the metacommunity species pool S. 

To reveal the origins of community types in the human microbiome, we decomposed 
the universal interaction matrix as 


A = NH o Gs, (4) 

which contains four components, (i) N G is the nominal interspecific interaction ma¬ 

trix where each element is sampled from a normal distribution with mean 0 and variance 
cr^, i.e. [N]ij ~ A/’(0, cr^). (ii) H G is a diagonal matrix that captures the overall 

interaction strength heterogeneity of different species. When studying the impact of inter¬ 
action strength heterogeneity the diagonal elements of H will be drawn from a power-law 
distribution with exponent —a, i.e. [H]^^ ~ 'P(«), which are subsequently normalized so 
that the mean of the diagonal elements is equal to 1. This is to ensure that the average 
interaction strength is bounded. For studies that do not involve interaction strength het¬ 
erogeneity H is simply the identity matrix, (iii) G G is the adjacency matrix of the 

underlying ecological network: [G]ij = 1 if species i is affected by the presence of species j 
and 0 otherwise. For details on the construction of G for different network topologies see 
Supplementary Text Sec. 3.2.2. Note that the Hadamard product (o) between H and G 
represents element-wise matrix multiplication, (iv) The last component s is simply a scal¬ 
ing factor between 0 and 1. Finally, we set [A]^^ = — 1. The presence of the scaling factor 
s and setting the diagonal elements of A to —1 are to ensure an asymptotic stability con¬ 
dition for the GLV dynamics (Supplementary Text Sec. 4.2, 4.3.3, and 4.5). The elements 
in the global growth rate vector r are taken from the uniform distribution, [r]^ ~ Z//(0,1). 
Details concerning the distribution A/", V and lA can be found in Supplementary Text Sec. 
3.1.1. 
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4 Origins of Community Types 

We first studied the role of interspecific interaction strength heterogeneity on the emergence 
of community types. In order to achieve this, we chose the complete graph topology, i.e. 
each species interacts with all other species. This eliminates any structural heterogeneity. 
The nominal interaction strengths were taken from a normal distribution A/”(0,1), the 
scaling component was set to s = 0.7, and the interaction strength heterogeneity was 
varied from low heterogeneity (a = 7) to a high level of heterogeneity (a = 1.01). Figure 1 
displays the distributions of the diagonal elements of the interaction heterogeneity matrix 
H at various heterogeneity levels. For each level of heterogeneity we constructed 500 LCs, 
each with 80 species randomly drawn from a metacommunity of 100 species. Figure lb 
illustrates the global interaction matrix A as a weighted network. With low heterogeneity 
all the link weights are of the same order of magnitude. As the heterogeneity increases 
fewer nodes contain highly weighted links, until there is only one node with highly weighted 
links when a = 1.01. These nodes with highly weighted links correspond to SISs. 

Figure Ic presents the results of Principle Coordinates Analysis (PCoA) of the steady 
states associated with the 500 different LCs as a function of a. For low interaction het¬ 
erogeneity (a = 7) the classical clustering measure. Silhouette Index, is less than 0.1, 
suggesting a lack of clustering in the data. As the heterogeneity increases the steady 
states can be seen to separate in the first two principle coordinate axes. At one point 
{a = 2.0) three clusters is the optimal number of clusters. Then as a continues to decrease 
the optimal number of clusters becomes two. The fact that there are three clusters when 
a = 2.0 is not special, as a different number of optimal clusters can be observed with 
different model parameters or different clustering measures (see Supplementary Text Sec. 
7.2) [52]. While the precise number of clusters is not important here, what is important is 
the fact that the degree of interaction strength heterogeneity controls the degree to which 
the clusters appear to be distinct. For low levels of interaction strength heterogeneity 
the clusters appear to be more like dense regions within a continuous gradient. As the 
heterogeneity increases, the clusters become more distinct. Indeed, having two clusters for 
a = 1.01 is to be expected, because one of the clusters is associated with all the LCs that 
contain the single SIS, and the other LCs that do not contain the single SIS constitute the 
other cluster. 

The overall trend observed in Figure Ic is unaffected if the complete graph is replaced 
by an Erdos-Renyi (ER) random graph, or if the total number of LCs is increased (Supple¬ 
mentary Figures SI and S2). The result is also generally unaffected by the specifics of the 
nominal distribution (Supplementary Text Sec. 7.2.1), the mean degree of the ER graph 
(Supplementary Text Sec. 7.2.2), or the number of species in the LCs (Supplementary 
Text Sec. 7.2.3). Of course, each LC can be invaded by other species that are currently 
absent. If this migration occurs relatively fast, then all LCs will converge to roughly the 
same species collection and the clustering will disappear. Hence in our modeling approach 
we have to assume that the migration occurs at a relatively slow time scale, and the time 
interval between species invasions is too long to disrupt the clustering. We also note that 
if heterogeneous interactions are placed at random in the network the clustering of steady 
states does not arise (Supplementary Eigure S3). Our results are also robust (in the control 
theoretical sense) to stochasticity and the migration of existing species [34] . Robustness to 
migration is illustrated in Supplementary Eigures S4 and S5, and robustness to stochastic 
disturbances is illustrated in Supplementary Eigures S6-S8 (see Supplementary Text Sec. 
4.4 for analytical robustness results). 

We can explain the above results as follows: for low interaction strength heterogeneity 
all of the matrices are very similar. In other words, despite containing different sets of 
species, all the LCs have very similar dynamics. Thus, clustering of steady states is not to 
be expected. As the heterogeneity of interaction strength increases, however, some of the 
LCs will have species that are associated with the highly weighted columns in A, i.e. the 
SISs. Eigure 2 presents a detailed analysis of the most abundant (dominating) species in 
each of the three clusters (community types) in Eigure Ic for a = 2 and a = 1.6, along with 
the abundances of the SISs within each cluster. It is clear that for different clusters their 
dominating species are different, consistent with the empirical finding that each enterotype 
is dominated by a different genus [3]. The SISs that are present in each cluster also vary. 
Eor instance with a = 1.6 all LCs in the blue cluster contain SISs number 23 and 81, and 
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none have species 60 or 51. For the orange cluster it is the opposite scenario. All of the 
LCs in the orange cluster contain SISs 60 and 51, and do not contain species 23 or 81. 
Most of the LCs in the yellow cluster contain SISs 23 and 51. Hence, each community 
type is well characterized by a unique combination of SISs. Note that none of the SISs 
are dominating species. These findings, along with the analysis in Box 1 , suggest that 
heterogeneity in interaction strengths or the presence of SISs leads to the clustering of 
steady states, i.e. the emergence of community types. 

We then studied the impact of structural heterogeneity on community types. Four 
different scenarios are illustrated in Figure 3: (a) a complete graph topology as in Figure 
1; (b) an ER random graph as in Supplementary Figure 1; (c) a power-law out-degree 
network; (d) a power-law out-degree network with no interaction strength heterogeneity. 
Figures 3a, 3b and 3c support the main result shown in Figure 1, i.e. increasing interac¬ 
tion strength heterogeneity leads to the emergence of distinct community types. Figure 
3d displays rather unexpected results as it suggests that structural heterogeneity alone 
does not lead to distinct community types. It is only with the inclusion of interaction 
strength heterogeneity that structurally heterogeneous microbial ecosystems can display 
strong clustering in their steady states as shown in Figure 3c. This result is rather sur¬ 
prising, because structural heterogeneity is observed in many real-world complex networks 
[1,6,58] and has been shown to affect many dynamical processes over complex networks 
[60,72,75]. 

Note that in the preparation of Figure 3 the steady state abundances were normalized 
to get relative abundances of the species and the Jensen-Shannon distance metric was 
used for clustering analysis [39]. The trends discussed above also hold when, instead of 
the Silhouette Index, the Variance Ratio Criterion is used as the clustering measure, or 
the Euclidean distance is used for clustering, or when absolute abundances are analyzed 
along with the Euclidean distance being used (Supplementary Eigures S9, SIO, and SI I). 
Supplementary Eigure SI I correlates to the analytical results in Box I , where absolute 
abundances and the Euclidean distance are implicitly used. 

Control of Community Types 

With the knowledge that each community type can be associated with a specific collection 
of SISs, we tested the hypothesis that a local community could be steered to a desired 
community type by controlling the combination of SISs only. Our results for three different 
scenarios are shown in Eigure 4a for a = 1.6. The local community that was controlled 
in each scenario is shown in magenta and is denoted LC*, which initially belongs to the 
blue cluster. Eor Scenario I, LC* had the SISs 23 and 81 removed, with species 60 and 
51 simultaneously introduced with random initial abundances drawn from Z4(0,1). Recall 
that species 60 and 51 are the SISs present in the orange cluster. This swap of SISs shifts 
LC* to a slightly different state (green dot) within the blue cluster. The GLV dynamics 
were then simulated and the trajectory goes from the blue cluster to the orange cluster. 
This result was independent of the initial condition of species 60 and 51 (Eig. 4b). This 
open-loop control of the community type by manipulating a set of SISs also works at lower 
levels of heterogeneity (Eig. 4c and 4d). Here we use the term open-loop to contrast 
closed-loop control where inputs are designed with feedback so as to continuously correct 
the system of interest. These findings imply that the SISs, despite their low abundances, 
can be used to effectively control a microbial community to a desired community type. 

In Scenario 2 we tested if the same result could be obtained by removing the six 
most abundant species from LC* and introducing the six most abundant species from the 
orange cluster at exactly the same abundance level as an arbitrary local community in 
the orange cluster. The state after this dominating species swap (red dot) starts close to 
the orange cluster, because the six most abundant species from a local community in that 
cluster were copied. The trajectory does not ultimately converge near the orange cluster, 
but goes toward the blue cluster instead. The trajectory, however, does not ultimately 
converge in the blue cluster because it does not contain any of the most abundant species 
present in the blue cluster. 

In scenario 3 we explored how the open-loop control methodology just presented could 
also be used to conceptually justify the success of EMT in treating patients with rCDI 
[98,105,106]. This scenario begins by removing 20 species from LC* (the top two SISs and 
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18 of the most abundant spaces) so as to emulate the effect of broad-spectrum antibiotics, 
resulting in an altered community (blue dot). Then the GLV dynamics were simulated and 
the local community converged to a new steady state (black dot), representing the GDI 
state. To emulate an oral capsule FMT 1% of the species abundances from an arbitrary 
LG in the orange cluster, i.e. the donor, was added to the GDI state, resulting in a slightly 
altered community (gray dot). The GLV dynamics were then simulated until the hnal 
steady state was reached (white dot). As expected the post-FMT steady state is in the 
orange cluster, the same cluster that is associated with the donor’s LG. Note that if during 
the FMT the SISs in the donor’s LG were not transplanted then the patient’s post-FMT 
steady state does not converge in the orange cluster (Supplementary Figure S12). 

The above results indicate that the presence of SISs simplihes the open-loop control 
design. However, the existence of community types is not a prerequisite for deploying 
this control methodology. The possibility for open-loop control of the human microbiome 
will likely be body site specific. Our work focused on the gut specifically because of the 
fact that this microbial community is very likely dominated by microbe-microbe and/or 
host-microbe interactions, rather than external disturbances. It is yet to be determined 
what factors drive the dynamics in other body sites. 

5 Discussion 

In this work we studied compositional shift as a function of species collection using a 
dynamic systems approach, aiming to offer a possible mechanism for the origins of com¬ 
munity types. We found that the presence of interaction strength heterogeneity or SISs is 
sufficient to explain the emergence of community types in the human microbiome, indepen¬ 
dent of the topology of the underlying ecological network. The presence of heterogeneity 
in the interspecific interaction strengths in natural communities has been well studied in 
macroecology [26,27,64,73]. Extensive studies are still required to explore this interesting 
direction in the human microbiome. While preliminary analysis is promising, all existing 
temporal metagenomic datasets are simply not sufficiently rich to infer the interspecific 
interaction strengths among all of the microbes present in and on our bodies [32] even 
at the genus level, let alone the species level. Recent studies have tried to overcome this 
issue by only investigating the interactions between the most abundant species [33]. Our 
results, however, suggest that SISs need not be the most abundant ones and can still play 
an important role in shaping the steady states of microbial ecosystems. Ignoring the lack of 
sufficient richness, system identification analysis with regularization and cross-validation 
[11,86] of the largest temporal metagenomic dataset to date [13] does not disprove the 
existence of SISs. To the contrary, it supports this assertion (see Supplementary Figure 
SI 3). Permutation of the time series however also results in the identification of interaction 
strength heterogeneity (see Supplementary Figures SI4 and SI5). Hence, the presence of 
SISs needs to be systematically studied with novel system identification methods and per¬ 
haps further validated with co-culture experiments [32]. For example, we could first use 
metabolic network models to predict levels of competition and complementarity among 
species [56], which could then be used as prior information to further improve system 
identification [2]. 

Note that our notion of SIS is fundamentally different from that of keystone species, 
which are typically understood as species that have a disproportionate deleterious effect 
(relative to its abundance) on the community upon their removal [74]. One can apply a 
brute-force leave-one-out strategy to evaluate the “degree of keystoneness” of any species 
in a given community network [7]. Even without any interaction strength heterogeneity, 
a given community may still have a few keystone species. The SISs defined in this work 
are those species that have very strong impacts (either positive or negative) on the species 
that they directly interact with. The presence of SISs requires the presence of interaction 
strength heterogeneity. We emphasize that an SIS is not necessarily a keystone species. 
In fact, without any special structure embedded in the interaction matrix (and hence the 
ecological network), there is no reason why the removal of any SIS would cause mass 
extinction. It does have a profound impact on the steady-state shift, which is exactly 
what we expected from our analytical results presented in Box I. 

Our findings also have important implications as we move forward with developing 
microbiome-based therapies, whether it be through drastic diet changes, EMT, drugs, or 
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even engineered microbes [53,55,65,85,89,103,104]. Indeed, our results suggest that a few 
strongly interacting microbes can determine the steady state landscape of the whole mi¬ 
crobial community. Therefore, it may be possible to control the microbiome efficiently by 
controlling the collection of SISs present in a patient’s gut. Finer control may be possible 
through the engineering of microbes. This will involve a detailed mechanistic understand¬ 
ing of the metabolic pathways associated with the microbes of interest. As discussed in 
Box 1 , given a new steady state of interest, the parameters 5, c, d, s could be chosen such 
that the new steady state is feasible and stable (Supplementary Text Sec. 4.3.1). Then, 
with the knowledge of the appropriate parameters 5, c, d, s it would be possible to introduce 
a known microbe with those characteristics or engineer one to have the desired properties. 
We emphasize that the stability and control of the microbial ecosystem must be studied 
at the macroscopic scale using a systems and control theoretic approach. This is similar 
to what is carried out in aerospace applications. The design of wings and control surfaces 
for an aircraft incorporate sophisticated fluid dynamic models. The control algorithms for 
planes however are often derived from simple linearized reduced order dynamic models 
where linear control techniques can be easily deployed [87]. Taken together, our results 
indicate that the origins and control of community types in the human microbiome can 
be explored analytically if we combine the tools of dynamic systems and control theory, 
opening new avenues to translational applications of the human microbiome. 
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Box 1. Steady State Shift in the Generalized Lotka-Volterra Model 


Impact of one new species: Consider the addition of one new species to the system described by 
(2), so that now we have m = n + 1 species, with the abundance of the z-th species denoted as 
Ziit), i = • , m. The new m-dimensional state vector z{t) = [z\{f), 2 : 2 (t), . . . , Zm{t)V evolves 


as 


z{f) = diag( 2 ;(f))(^ + Fz{f)), 


(Bl) 


where 



and 




Note here that A and r are as defined in (2) and we have only introduced four new elements 
s,b,c,d. The scalar element s represents the growth rate of the additional species Zm- The 
n-dimensional vector b represents the impact that species Zm has on the first n species zi:n, 
which also correspond to the n species in the original state vector x. The scalar element d is the 
Verhulst term for the new species Zm- Finally, the n-dimensional vector c represents the effect 
that zi:ri have on the dynamics of the m-th state Zm- 


Given the dynamics in (2) and (Bl), with the assumption that the matrix A is full rank, we can 
show that the difference between the original steady state x* and the shifted steady state of the 
same n species, denoted as z^.^, satisfies the following equality 

(B2) 

where z^ is the steady state value of the newly added species. Given that A is the same in (2) 
and (Bl), the shift in the steady state of the original n species is bilinear in terms of the vector b 
and the steady state value of the newly added m-th species z^. 

different initial original steady state 



Open-loop Control of Steady State: If there exists a diagonal matrix P > 0 such that 

A^P + PA < 0 then for any 2 ;* there exists 6, c, d, s such that the original steady state x* of (2) 
can be steered to 2 ;* and furthermore this 2 ;* can be made to be uniformly asymptotically stable 
(Supplementary Text Sec. 4.3.1, Theorem 8). 

Impact of multiple non-common species: Gonsider two systems with different collections of species. 
Let be the steady state of system 1 and z* be the steady state of system 2. Assuming that the 
systems share n species, then we can just apply the results from (B2) recursively to calculate the 
difference between the abundances of the common species as 

- n*. = (e - E ^^4] 

\ieM ieM / 

where M and M are the indices of the non-common species of system 1 and system 2 respectively. 
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a 



Principal Coordinate 1 Principal Coordinate 1 Principal Coordinate 1 Principal Coordinate 1 Principal Coordinate 1 


Figure 1: Impact of interaction strength heterogeneity on the distinctness of 
community types. A total of q = 500 local communities, each with p = 80 species 
randomly drawn from a universal pool of n = 100 species. The nominal components 
were drawn from Af{0, 1), the interaction heterogeneity matrix elements were taken 
from V{a) and a is varied with the set of values {7, 3, 2, 1.6, 1.01} for each 
column in the figure. The topology component G has all elements equal to 1, giving 
a complete graph. The scaling factor was set at s = 0.07. (a) Histogram of the 
diagonal elements of the heterogeneity matrix H. (b) Visualization of the universal 
interaction matrix A as a weighted adjacency matrix of a digraph, (c) Principle 
coordinate analyses of the normalized steady state for each local community using 
the Jensen-Shannon distance. The Silhouette Index and optimal number of clusters 
are denoted. Further details can be found in the Methods Section. 
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Figure 2: Comparison of dominating species to SISs in different community 
types (clusters). The relative abundances of the six most abundant species from 
each of the three clusters in Figure Ic for a = 2 and a = 1.6 are compared to that 
of the four species with the largest interaction strengths (60, 23, 81, and 51). 
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Figure 3: Impact of network structure on the distinctness of eommunity types. 
For each type of network structure 10 different Universal Triples (S, A, r) with n = 
100 species and q = 500 local communities of size p = 80 were generated with results 
shown in the lighter color and averaged results shown in bold, (a) Complete graph. 
Same study as in Figure 2 with a G [5, 1). (b) Erdos-Renyi network (digraph) 

[N]ij ~ A/’(0, 1), [H.]ii ~ Vioi) where a G (1,5], Probability [G]ij = 1 is 0.1, i.e. a 
mean in(out)-degree of 10, and scaling factor s = l/\/T6. (c) Power-law out-degree 
network [Nj^j ~ A/’(0, 1), [Hj^^ ~ V{o!.), G is the adjacency matrix for a digraph 
with out-degree having a power-law distribution V{a). The high-degree nodes have 
the largest interaction scaling, (d) Power-law out-degree network, no interaetions 
strength heterogeneity [Nj^^ ~ A/'(0, 1), H is the identity matrix, G is the adjacency 
matrix for a digraph with out-degree having a distribution V(a). Further details 
can be found in the Methods Section. 
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Original community LC* 
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interacting species swap 
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antibiotics 
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Figure 4: Open-loop control of the human microhiome. (a) Background of clus¬ 
tering analysis for a = 1.6 from Figure 1, but with Euclidean distance used so that 
a projection matrix could be found to show the trajectories in the 2D principle co¬ 
ordinate plane (Supplementary Text Sec. 5.6). We aim to steer a local community 
(denoted as LC*, shown in magenta) in the blue cluster to the orange cluster. Three 
different scenarios are presented, per the three numbers above the arrows. Scenario 
1: SISs swap. The SISs (23 and 81) of LC* were replaced by the SISs present in the 
orange cluster (60 and 51). The initial abundances of species 60 and 51 were drawn 
from U{0, 1), resulting in the green dot, and the GLV dynamics were simulated. 
Scenario 2: dominating species swap. The six most abundant species in LC* were 
removed and replaced by the six most abundant species from a local community 
in the orange cluster, with the initial condition after the switch of species shown 
as the red dot, and the dynamics were simulated until steady state was reached. 
Scenario 3: Fecal Microbiota Transplantation (FMT). The two SISs and 18 of the 
most abundant species (for a total of 20) were removed from LC* with the initial 
condition shown in blue (post-antibiotic state). Then the GLV dynamics were sim¬ 
ulated (gray line) and the system converged to the black dot (GDI state). Then 1% 
of the steady abundances from an arbitrary LG in the orange cluster were added to 
the GDI state (gray dot, emulating oral capsule FMT) and the dynamics were then 
simulated until steady state was reached, (b) The SISs swap process was repeated 
ten times, each time the initial abundances of species 60 and 51 were randomly 
drawn from U(0, 1). Nine of the simulations are shown in black and the simulation 
that pertains to Figure 4a is shown in maroon, (c) The same analysis as for Figure 
4a, in terms of SISs swap, but for a = 2. (d) The same analysis as for Figure 4a, 
in terms of SISs swap, but for a = 3. 
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Chapter 3 


Materials and Methods 


The methods section begins with a toy example to illustrate the construction of the uni¬ 
versal interaction matrix A = NH o Gs in (4), where 


steps: (i) 
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Given that H is diagonal, it scales the columns of N. If one thinks of A as the adjacency 
matrix of a digraph, then H scales all of the edges leaving a node. Thus one can con¬ 
sider H as controlling the interaction strength heterogeneity of A. Given the Hadamard 
product between H and G, the off-diagonal elements of G that are zero will result in the 
corresponding off-diagonal elements of A being zero as well. 

In the first study (Figure 1), to explore the impact of interaction heterogeneity on 
steady state shift, we varied the exponent —a of the power-law distribution of [H]ii to 
generate five different universal interaction matrices A of dimension 100 x 100. For each 
universal interaction matrix A, the nominal component N consists of independent and 
identically distributed elements sampled from a normal distribution A/'(0,1). The topology 
for this study was a complete graph and thus all the elements in G are equal to 1. The 
heterogeneity element H is constructed in two steps. First, five different vectors h{a) G 
are constructed where each element is sampled from a power-law distribution V{a) 
for a G {7,3,1.6,1.2,1.01}. Then, each of the h{a) is normalized to have a mean of 1, 
h = h/meaniji). Finally the heterogeneity matrix is defined as H = diag(/i). For this 
study s = 0.07, ensuring uniform asymptotic stability for the case of low heterogeneity 
(see Supplementary Text Theorem 17). The final step in the construction of A is to set 
the diagonal elements to —1. 

For each a the following simulation steps were taken. There are a total of 100 species, 
S = {1, 2, ..., 100}, in the metacommunity, and each of the 500 local communities con¬ 
tains 80 species, randomly chosen from S. The MATLAB command used to perform this 
step is randperm. The initial condition for each of the 500 local communities, x^^^(O), 
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were sampled from Vl{0, 1). The dynamics were then simulated for 10 seconds using the 
MATLAB command ode45. If any of the 500 simulations crashed due to instability or if 
the norm of the terminal discrete time derivative was greater than 0.01 then that local 
community was excluded from the rest of the study. Those simulations that finished with¬ 
out crashing and with small terminal discrete time derivative were deemed steady. Less 
than 1% of simulations were deemed unstable in the preparation of Figure 1. 

The networks presented in the second row of Figure 2 were constructed by considering A 
as the weighted adjacency matrix of the network. Note that arrows showing directionality 
and self loops were suppressed. The links were color coded in proportion to the absolute 
value of the entries in A. 

For the last row of Figure 1 a clustering analysis was performed. For each a the steady 
state abundances of the 500 local communities were normalized so that we have 500 syn¬ 
thetic microbial samples. Then /c-medoids clustering was performed for /c G {1, 2, ..., 10} 
using the Jensen-Shannon distance metric (Supplementary Text Sec. 5.1). Silhouette anal¬ 
ysis was performed to determine the optimal number of clusters and the clustering results 
were illustrated in the 2-dimensional principle coordinates plot. For Supplementary Figure 
SI the same steps as for the preparation of Figure 1 were performed, but with G repre¬ 
senting the adjacency matrix of an Erdos-Renyi digraph with mean degree of 20 (mean 
in-degree of 10 and mean out-degree of 10) and s = 1/a/TO- Details on the construction 
of an Erdos-Renyi digraph can be found in Supplemental Information Section 3.2.1. Eor 
Supplementary Eigure S2 the same steps as above were performed in Eigure 1 but with 
p = 5, 000 local communities. 

Eigure 3 is a macroscopic analysis of how network structure plays a role in the steady 
state shift with values of a G (1,5]. Eor each topology ten different universal matrices A 
were generated. Eigure 3(a) shows the results of a complete graph and for each of the ten 
universal A the same steps as in the preparation of Eigure 1 were carried out. Eigure 3(b) 
shoes the result of an Erdds-Renyi random digraph topology and for each of the ten A 
matrices the same steps as in the preparation of Supplementary Eigure S2 were carried 
out. Eigure 3(c) shows results for networks with a power-law out-degree distribution with a 
mean out-degree of 10, where the out-degree sequence uses the same h in the construction 
of H. More information on the construction of G for a power-law out-degree network can 
be found in Supplementary Text Sec. 3.2.2. Eigure 3(d) shows results for networks with a 
power-law out-degree distribution with mean out-degree of 10 and there is no interaction 
strength heterogeneity, i.e. H is the identity matrix. Eor this study the Silhouette Index 
was constructed from normalized steady state data using the Jensen-Shannon distance. 
Supplementary Eigure S9 is the same as Eigure 3, but instead of the Silhouette Index, 
the variance ratio criterion is used with the Jensen-Shannon distance, from normalized 
steady state abundance (Supplementary Text Sec. 5.4). In Supplementary Eigure SIO the 
Silhouette Index is determined from the Euclidean distance with normalized steady state 
abundance. Einally, in Supplementary Eigure Sll the Silhouette Index is determined by 
the Euclidean norm with the absolute steady state abundance. 

Eigure 4 contains a PCoA analysis of the results from Eigure 1, but with the Euclidean 
distance being used instead of the Jensen-Shannon distance, making PCoA equivalent to 
principle component analysis. This enables us to project the open-loop control trajectories 
into the principle coordinates (Supplementary Text Sec 5.6). This procedure was also used 
in the preparation of Supplementary Eigure S12 . 

Supplementary Eigures S13 to S15 contain system identification analyses for temporal 
gut microbiome data of two subjects [13]. The data is publicly available from the metage¬ 
nomics analysis server MG-RAST:4457768.3-4459735.3 and can also be accessed (as we 
did) from Qiita (http://qiita.ucsd.edu) under study ID 550. The processed data was 
downloaded as biom file “67_otu_table.biom” (2014-11-17 13:18:50.591389). The Opera¬ 
tional Taxonomic Units (OTUs) were then grouped from the genus level and up, depending 
on the availability of known classifications for OTUs, and converted to a txt file using Mac- 
QIIME version 1.9.0-20140227 with the command suinmarize_taxa.py with the options 
-L 6 -a true. Data was collected over 445 days with 336 fecal samples from Subject A 
and 131 fecal samples from Subject B. Details on the system identihcation algorithm are 
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now given. The dynamics in (2) can be approximated in discrete time as [86] 

n 

ei{k) + \og{xi{tk+i)) - \og{xi{tk)) = n + '^aijXj{tk) (Ml) 

j=i 

for z = 1, 2,..., n where /c = 1, 2,..., iV — 1 is the sample index, N is the total number of 
samples, tk is the time stamp of sample /c, and e is an error term that arises because of 
the assumption that x{t) is constant over each interval t G Equation (Ml) can 

be rewritten in terms of a regressor vector 

0(/c) = [l,Xi(tfc), X2{tk), Xn{tk)]^, 

the parameter vector Oi = [r^, an, ai 2 , ..., cizn] and the log difference yi{k) = log {xi(tk-\-i)) — 
log {xi(tk)) as 

ei{k) + yi{k) = Oi(l){k). 

The identification problem can then be defined as finding the parameter matrix estimate 
© = pi, Oln ■ ■ , of the true parameter matrix 0 = [^i, ? •'' ? ^1]^- Letting 

2/(^) = [2/i(^). ...,yn{k)Y 

be the log difference vector for all species and Y — [z/(l), ?/(2), ..., y{N — 1)] be the log 
difference matrix the system identification problem can be compactly presented as 

min||y-0$||| + A||e||| 

© 

where $ = [0(1)? 0(2), ..., 0(iV — 1)] is the regressor matrix, \\-\\f denotes the Frobenius 
norm, A > 0 is the Tikhonov regularization term [96]. The minimal solution to the above 
problem can be given directly as 

argmin (||r - 0-I>|||. + A||0|||.) = FY^{YY^ + A/)"^ 
where I is the identity matrix. 

Next we discuss how missing data, zero reads, and A were chosen. The difference 
equation in (Ml) only uses sample data over two consecutive time samples. Therefore, in 
the construction of Y and $ we only include samples that for which there is data from the 
next day as well. Also, given that the logarithms are used, when a sample has zero reads for 
a given taxa, a read value of one is inserted. Then relative abundances are computed before 
the logarithm is taken. Finally we discuss how the regularization parameter is chosen. 
For Supplementary Figures S13 and S14 the following cross-validation is performed. For 
Subjects A and B two-thirds of data was used for training and one-third for testing. More 
precisely, for each A two-thirds of the data from Subject A and two-thirds of the data 
from Subject B were used to identify their corresponding dynamical constants. Then the 
combined error from the two test sets was used to find the optimal A. The regularization 
value used in Supplementary Figure S15 is simply the same regularization value used in 
Supplementary Figure S13 . 
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Supplementary Figures 
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Supplementary Figure SI: Impact of interaction strength heterogeneity on 
the distinctness of eommunity types. Same as Figure 1 but with the topology 
component G chosen to be an Erdos-Renyi digraph with a link probability of 0.1 
and the scaling factor was set at s = 1/VTo. 
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Principal Coordinate 1 Principal Coordinate 1 Principal Coordinate 1 Principal Coordinate 1 Principal Coordinate 1 


Supplementary Figure S2: Impact of interaction strength heterogeneity on the 
distinctness of eommunity types. Same as Figure 1 but with p =5,000 local com¬ 
munities. Note that it is rather counter-intuitive that for a = 1.01 the Silhouette 
Index suggests that there are two clusters, while PCoA suggests three clusters. We 
emphasize that as a typical ordination method, the PCoA just produces a spatial 
representation of the entities in the dataset, rather than the actual determination 
of cluster membership [44,67]. Note that as compared to Figure 1, because there 
are more samples in this figure, the distinctness of the clusters when a = 2 has 
shifted to more of a continuous gradient as apposed to distinct clusters. 
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Supplementary Figure S3: Impact of interaction heterogeneity disbursed ran¬ 
domly throughout the network. The set up is the same as that of Figure 1 but 
instead of H being a diagonal matrix, it is a full matrix, so that individual interac¬ 
tions are scaled randomly from a power-law distribution. 
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Supplementary Figure S4: Impact of low levels of migration. Same as Figure 
1 but with a new term A(t) G [0, added to the dynamics so that now x = 
A + diag(tc)(r + Ax). In this example \i ~ ^/(0,0.1). The disturbance is sampled 
every 0.01 seconds and held constant until the next sample is taken.. 
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Supplementary Figure S5: Impact of moderate levels of migration. Same as 
Figure 1 but with a new term A(t) G [0, 1]^ added to the dynamics so that now 
X = A(t) + diag(tc)(r + Ax). In this example Xi ~ ^/(O, 1). The disturbance is 
sampled every 0.01 seconds and held constant until the next sample is taken. 
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Supplementary Figure S6: Impact of small stochastic disturbance. Same as 
Figure 1 but with stochastic Ito dynamics dx = diag(tc)(r dt + Ax dt + cdw) where 
It; is a n-dimensional Brownian motion and c represents the stochastic disturbance 
strength. Dynamics were simulated with a discrete time step of 0.01 seconds and 
c = 0.1. 
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Supplementary Figure S7: Impact of moderate stochastic disturbance. Same as 
Figure 1 but with stochastic Ito dynamics dx = diag(tc)(r dt + Ax dt + cdw) where 
It; is a n-dimensional Brownian motion and c represents the stochastic disturbance 
strength. Dynamics were simulated with a discrete time step of 0.01 seconds and 
c = 0.5. 
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Supplementary Figure S8: Impact of large stochastic disturbance. Same as 
Figure 1 but with stochastic Ito dynamics dx = diag(tc)(r dt + Ax dt + cdw) where 
It; is a n-dimensional Brownian motion and c represents the stochastic disturbance 
strength. Dynamics were simulated with a discrete time step of 0.01 seconds and 
c = 1. 
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Supplementary Figure S9: Impact of network structure on the distinctness 
of eommunity types. The same as Figure 3 with the Varianee Ratio Criterion 
(VRC) used as apposed to the Silhouette Index for the clustering measure. See 
Supplemental Information §5.4 for details on the VRC. 
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Supplementary Figure SIO: Impact of network structure on the distinctness of 
community types. The same as Figure 3 with the Euclidean distance metric used 
instead of the Jensen-Shannon distance metric. 
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Supplementary Figure Sll: Impact of network structure on the distinctness 
of eommunity types. The same as Figure 3 with the Euclidean distance metric 
used instead of the Jensen-Shannon distance metric and absolute abundance used 
instead of relative abundace. 
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Supplementary Figure S12: Unsuccessful Fecal Microbiota Transplantation. 
Similar to Scenario 3 shown in Figure 4a, but during the FMT, the SISs (60 and 
51) of the donor’s local community in the orange cluster were not transplanted to 
the CDI state (black dot). This FMT resulted in a slightly altered community (gray 
dot) and the system eventually evolved to a steady state (white dot) which is not 
in the orange cluster. Hence the FMT failed. 
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Supplementary Figure S13: System Identification, Tikhonov Regularization 
A = 0.0423. System identification was performed on the stool samples from the 
longitudinal data in [13] for two subjects as described in the Supplementary Meth¬ 
ods where A was determined by cross-validation, (a) Visualization of microbial taxa 
in terms of relative abundances versus day sample was taken, (b) Heat map of the 
interaction matrix for top 100 SISs. (c) Histogram of Standard Deviation (SD) of 
the columns of the interaction matrix, (d) List of top ten SISs in descending inter¬ 
action strength (defined by the SD of each column in the interaction matrix) with 
relative abundances over all samples shown as a box plot. The banded structure 
shown in the heat map supports the assertion that SISs do exist in the gut micro- 
biome. However this banded structure is also seen when the dates of the sample 
collections are permuted, see Supplementary Figure S14 and S15. 
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Supplementary Figure S14: System Identification, Day Swap, Tikhonov Reg¬ 
ularization A = 0.0057. System identification was performed on the stool samples 
from the longitudinal data in [13], but with the collection dates permuted, A was 
determined by cross-validation on the permuted data, (a) Visualization of microbial 
taxa in terms of relative abundances versus day sample was taken (not permuted 
samples), (b) Heat map of the interaction matrix for top 100 SISs. (c) Histogram 
of Standard Deviation (SD) of the columns of the interaction matrix, (d) List of 
top ten SISs in descending interaction strength (defined by the SD of each column 
in the interaction matrix) with relative abundances over all samples shown as a box 
plot. Even though the sample days have been permuted the banded structure still 
persists. 
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Supplementary Figure S15: System Identification, Day Swap, Tikhonov Reg¬ 
ularization A = 0.0423. System identification was performed on the stool samples 
from the longitudinal data in [13], but with the collection dates permuted, A was 
selected to be the same as in Supplementary Figure S13. (a) Visualization of mi¬ 
crobial taxa in terms of relative abundances versus day sample was taken (not 
permuted samples), (b) Heat map of the interaction matrix for top 100 SISs. (c) 
Histogram of Standard Deviation (SD) of the columns of the interaction matrix, 
(d) List of top ten SISs in descending interaction strength (defined by the SD of 
each column in the interaction matrix) with relative abundances over all samples 
shown as a box plot. For the permuted data when A is larger than the optimal value 
from the cross-validation the identification method biases towards making the most 
abundant species also the SISs. 
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Chapter 5 


Supplementary Text 


1 Introduction 

The purpose of this supplement is to give a self-contained, rigorous account of the topics 
needed to fully understand the main text. The supplement is organized as follows. Section 

2 contains some remarks on notation. Section 3 covers the basics of random variables, gives 
details on how one constructs a well defined random variable from a power-law distribution, 
and covers some basic results in random matrix theory, Wigner’s circle and semi-circle 
laws. Section 4 contains details on the definition of asymptotic stability, how this relates 
to the generalized Lotka-Volterra dynamics, some resent results on diagonal stability are 
discussed, and new results regarding the stochastic stability of Lotka-Volterra dynamics 
is discussed. Section 5 gives the details of the multi-dimensional scaling analysis and 
clustering algorithms employed in this study. Section 6 gives the details of our modeling 
approach. Section 7 contains simulation results that are intended to complement the 
hgures in the main text and the supplementary hgures. 

Much of the foundational material is well known in each of their respective areas. The 
diagonal stability result for random matrices has never been discussed in the context of 
Lotka-Volterra dynamics. The discussions regarding stable steady state shift for Lotka- 
Volterra dynamics are the first of their kind as well. 

Those familiar with probability theory need only read §3.3 for a refresher on Wigner 
random matrices and the rest of §3 can be skipped. Those familiar with stability need 
only visit §4.3 and §4.5. That being said, even the seasoned stability theorist may find 
new things in the stochastic subsection of §4.4. Section 5 can be skipped for those familiar 
with the commonly used tools in clustering analysis. 

2 Notation 

Throughout we denote the real numbers as R = (—oo, oc), and the complex numbers as C. 
A number zGCifz = x + zy where x, y G R and i = The positive real numbers are 

denoted as R>o = (0, oo) and the non-negative reals as R>o = [0,oo). The n-dimensional 
reals are denoted as R’^, R^o denotes the n-dimensional space of positive vectors which 
we will refer to as the positive orthant^ and R>o as the non-negative orthant. An m x n 
matrix A of values in the reals is denoted as A G The element of A in the z-th row 

and the j-ih column will often be denoted with lower case letters as Oij = [A]ij. If there 
is an issue with formatting and it is not clear that both elements i and j are subscripts in 
the previous notation then it is equivalent to denote Oij as aij. The dual notation for an 
n X n matrix constructed from elements aij is denoted as A = {aij)i<ij<n- 

The superscript (•)^ is used to denote transpose. The components of an n dimensional 
vector X are defined as follows x — [xi^ X 2 , ..., Xn]^■ We will often make use of the 
following nonstandard subscript notation, Xj:k to denote the vector obtained from taking 
the j-th element to the k-th element of As an example, consider y = [1, 3, 5]^, then 
2 / 2:3 = [3, 5]^. The Euclidean norm of a vector x G R’^ is defined as ||x|| = 

When applied to a matrix A G the norm is an induced norm ||A|| = sup|| 3 ,||^^ ||Ax||. 

^This notation was motivated by the index notation used in Matlab. 
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For a square symmetric matrix P — P^ the inequality < (<) is used as P < 

0 (P < 0) if and only if Px < 0 {x^Px < 0) for all x 

The notation for union and intersection of sets is U and H respectively. For a set A C U 
the complement of A is defined as = {y G P \y ^ A}. The set minus notation is defined 
as A \ B = An Let A be a set, then |^| denotes the cardinality of that set. If for 
instance A = {1, 2,4}, then \A\ = 3. 

Given that x will primarily be defined as the state variable, we will use Y when 
discussing a generic random variable. The probability distribution for Y is denoted as 
yy for which we will generically denote its probability density function as f{y) where 
fly = f f{y)dy. When given a probability density function /(y), the notation Y ~ f(y) 
reads as the random variable Y is drawn from the probability density function /. Also, 
if we generically denote the standard normal distribution as A/'(0,1), then with a slight 
abuse of notation we can write Y ~ ^"(0,1) to denote that Y is drawn from the standard 
normal distribution. The notation A = T is used to denote when two random variables 
are drawn from the same distribution. 

We will also make use of the following asymptotic notation, /(n) = 0{g{n)) if there 
exists a C independent of n such that ||/(n)|| < C||^(n)|| for n sufficiently large. Another 
form of asymptotic notation that will be borrowed is the following, f{n) = o{g{n)) if there 
exists a c{n) > 0 where limn^oo c(n) = 0 and ||/(n)|| < c(n)||^(n)||. 

Results from dynamics and control, to probability and random matrix theory will be 
called upon in this work. If the notation of a particular section seems to overlap, it 
is assumed that the convention from the field of origin overrides. For instance, when 
discussing dynamics a capital letter will denote a matrix, however capital letters are used 
in probability theory when denoting generic random variables. 

3 Random Variables and Random Matrices 
3.1 Primer on Random Variables 

We will use the following three distributions in constructing the random interaction ma¬ 
trices for the microbial communities: the uniform distribution taking values in the interval 
[0,1], the normal (Gaussian) distribution of mean 0 and variance and a power-law dis¬ 
tribution with minimum value 1 and exponent —a. One can define a random variable just 
with a variable space and ignore the underlying sample space [92, §1.1.2]. We however 
introduce the full probabilistic machinery into the discussion. Formality in this section is 
two fold. First, to give the reader confidence that the random variables generated from the 
power-law distribution are indeed well defined random variables [15]. Second, the Wigner 
circle law and semi-circle law can not be defined without this machinery, and these two 
laws are fundamental to understanding stability results that are presented in the next 
section. The following is a nearly verbatim presentation of probability spaces following 
[91,92]. 

Let Q be a sample space. When given a measure and a cr-algebra, Q becomes a 
probability space Q = (Q,A, P), where A is a cr-algebra of subsets of Q and P is a 
probability measure. The probability measure along with the sample space satisfy the 
following equality P(f^) = 1. Events E are then taken from the cr-algebra and will have 
a probability of occurring, i.e. E E T and E i— P(E) where P{E) G [0,1]. A random 
variable Y takes values in a measurable space R = (R, IZ) where 7^ is a cr-algebra of subsets 
of R. We will also refer to R as the variable space. The formal definition of a random 
variable is a map Y \ Pi ^ R where Y is measurable. When one asks for the probability 
that event T is in S' G 7^, we are interested in the probability of event E = Y~^{S) 
occurring. From the definition of our probability space Q, this is simply P(Y~^{S)). Note 
that this is equivalent to P({ct; G Pi : Y{u) G S}) which we can unambiguously denote in 
shorthand as P(y G S) [92, §1.1.2]. This notation makes no reference to the original 
sample space and thus can be used without ambiguity when discussing the probability of 
an event occurring. All of our probability spaces and variable spaces will be subsets of R 
or C. Thus we will implicitly use the Borel cr-algebra. Gonsequently, the specific cr-algebra 
of interest will no longer be denoted. 

We now give rigorous definitions for the probability measure gy of a random variable 
Y ^ R— (R,7V) and the corresponding probability density function / associated with gy. 
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The distribution of Y is defined as 


fiY{S)^R(Y eS). 

Given the above definition we define the probability density function as 

fiY{S)= [ f{y)dy. 

Js 

The cumulative distribution function is defined as 

Friy) ^P{Y < y) = 


The expected value of a random variable Y ^ f{y) (taking values from the probability 
density function /) is defined as 


ET = 


= [ ydviy) 

Jr 

/ yf{y)dy 

Jr 

[ PiY > X)dX. 

Jr 


(Tl) 


Sometimes parentheses are used, i.e. E(y), for clarity. The variance of a distribution is 
defined as 

Var(y) = E |y - E(y)|^ . (T2) 

Definition 1 ([Almost Surely). An event E occurs almost surely if P(^) = 1. 


3.1.1 Distributions of interest 

The uniform distribution generates random variables on the measure space R = [0,1] with 
the following probability density function 


u{y) 


1 if 2/ € [0,1] 
0 otherwise . 


A pseudorandom variable in [0,1] can be generated uniformly using the MATLAB com¬ 
mand random( ^ unif ^ ,0,1)- We will use the shorthand notation U(0, 1) to denote a uni¬ 
form distribution taking values in [0,1]. 

The normal distribution with mean 0 and standard deviation a generates random 
variables on the measure space R and satisfies the well known probability density function 


n(y) 



A pseudorandom variable with mean 0 and variance can be generated using the MAT- 
LAB command random(^norm^ ,0,sigma). We will use the shorthand notation 
to denote a normal distribution with mean 0 and standard deviation a. A log-normal 
distributed random variable is simply one in which Y = where X ~ A/'(0, a^). 

The generic power-law used in this work generates random variables on the measure 
space [1, oo) with the following probability density function 


p{y) = (a - 1)2/ 


(T3) 


where a > 1 [15, Equation (2.2)]. We will use the following short hand notation V{a) to 
denote a power-law distribution with exponent —a. 

Later we will generate power-law distributions from uniform distributions, and thus we 
have the following simple result. Let U he a, uniform random variable from the set [0,1), 
then we can generate a random variable P with a power-law distribution in [1, oo) using 
the following measurable monotonic function r : [0,1) ^ [1, oo). 


r{U) = (1-t7)T^ 


(T4) 
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Supplementary Text Figure Tl: Illustration of the distributions for the three 
types of probability functions used in this work, from left to right: uniform, normal, 
and power-law. 


Now, letting P = r{U) we can see that P indeed satisfies all of the requirements to be 
a well defined random variable. P is a measurable function defined from a probability 
space (the variable space of U with a probability measure P) to an event space by a 
measurable function r. This mapping for the generation of a random variable P with 
power law distribution is illustrated with the following commutative diagram. 


Q 



Random variables satisfying a power-law probability density function can be generated in 
MATLAB using the following line (l-random(^unif ^,0,1))"(1/(1-alpha)). 

The nice trick in (T4) is introduced in the literature without proof, and is simply 
denoted as following from the fundamental transformation law of probabilities [77, §7.3]. 
Given the rigorous introduction of probability theory at the beginning of this section, we 
can however directly derive this result by analyzing the probability that P > y under the 
assumption that there exists a monotonically increasing bijection r : [0,1) ^ [1, oo), which 
yields 


POO 

/ p(x)dx — P(P > y) 

J y 

= P(r(P) > y) 

= P{U>r-\y)) 


(T5) 




'-(y) 


u{x)dx. 


The transition from line 2 to line 3 in the above equality follows from the fact that r is 
monotonically increasing and thus when the inverse is taken the sign of the inequality is 
preserved. Integrating the above equality it follows that 


Substituting z = r ^(y) and solving for y in the above equality it follows that 


y = 


(1 — z) 


which proves the relation in (T4). 


3.2 Random Network Models 

Formally a digraph Q is defined by the double (V, 8) where V = {1, 2,..., n} is the vertex 
set and the directed edges are defined by the ordered pairs (z,j)G^CVxV. An element 
(i, j) G 8 if and only if there is a directed edge from vertex i to vertex j. We are primarily 
interested in the adjacency matrix A of the digraph which is defined as 

[Ai\ij — 


1 if (j, i) e 8 
0 otherwise 
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Power-law a = 7 Power-law a = 2 Power-law a = 1.2 





y y y 

Supplementary Text Figure T2: Illustration of the power-law distribution from 
low heterogeneity to high heterogeneity. 


When [A]ij = 1 for all i and j the digraph is said to be complete. Note that we allow for 
self-loops in our construction. Two other digraph topologies to be discussed shortly are 
the Erdds-Renyi (Gilbert) random digraph and the power-law degree digraph. 


3.2.1 Erdos-Renyi (Gilbert) Digraph 


An Erdds-Renyi (Gilbert) digraph is a digraph Q(n,p) of n nodes, where the probability 
of a directed edge from node i to node j is p for any z, j G V. The adjacency matrix 
for this model is constructed as follows. Let G G [0,be a matrix with elements 
independently sampled from the uniform distribution between 0 and 1, [G]ij ~ Zi(0,l). 
Then let A be defined as follows 


[A\ij 


1 if [G]ij < p 
0 otherwise 


If one is interested in defining an Erdds-Renyi model for a 100 node digraph with an 
expected mean in-degree (or out-degree) of 10, then simply set p = 10/100 in this con¬ 
struction. 

While the above model is often credited to Erdds and Renyi [28,29], it was actually 
first presented by Gilbert in [36]. We will follow convention however and simply refer to 
this random network model as the Erdds-Renyi (ER) model. The adjacency matrix for this 
model can be generated in MATLAB using the following boolean expression rand(n,n)<p. 


3.2.2 Power-law Out-degree Digraph 


In this section we outline how one can generate the adjacency matrix for a power-law 
out-degree digraph. Let h = [hi, ^ 2 , ..., hn]^ be the column vector of out-degrees for 
nodes {1, 2, ..., n}. If one is interested in having a digraph with a power-law out-degree 
of exponent —a with the mean out-degree approximately d, then setting 


hi = min 


d 


lh]i_ 

mean(h) 



(T6) 


is sufficient, where [h]i ~ V{a), i — 1,2,... ,n. Note that [•] is the ceiling operator. We 
also note that this will not guarantee that the mean out-degree is d for any hnite sized 
digraph. Einally the adjacency matrix is constructed by selecting hi random elements in 
column z of ^ and setting them to 1. 

This method of generating power-law degree distributions leaves much to be desired. 
It is not based on any known theory for power-law degree graphs [9,58]. So as to be able 
to compare Erdds-Renyi digraphs with the above power-law digraphs we have normalized 
by degree, which is not a common practice in the literature either. We note also that few 
authors have rigorously analyzed power-law digraphs with the exception of [9, §11] and 
[ 8 ], 


3.3 Spectrum of Random Matrices 

Two classic results from Wigner [99-101] are now discussed. Eirst we define the Empirical 
Spectral Distribution (ESD) of an n x n real valued matrix A as /za : C ^ N 

fiAiz) = — \{1 < i < n : ReXi{A) < Rez, ImAi(A) <lmz}\. 
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Recall that when applied to a finite set j-j denotes the cardinality of that set. The BSD 
simply counts the number of eigenvalues of A within radius jzj of the origin. A weaker 
version of [93, Theorem 1.10] is now stated. 

Theorem 1. Let Bn he a real n x n matrix whose elements are independently and identi- 
eally distributed random variables with mean 0 and varianee 1. Then it follows that Yn 
eonverges almost surely to the uniform disk in the eomplex plane, l|z|<i; with probability 
1. Let R denote the variable spaee for the i.i.d. random variables, then 


P [ limsup 

\ n^oo 


/ h{z)dnB^/^{z) - / h{z)li^i<i dz 
Jc Jc 


<e =1 


for all e > 0 and every bounded eontinuous funetion h : R ^ C. 


The semi-circle distribution is defined as /asc{S) = psc{y)dy where the probability density 
function is defined as 


Psc{y) 


\y\<2 

fo, \y\ > 2 


We now state a more conservative version of [90, Theorem 5]. 


Theorem 2. Let Mn be a symmetrie nxn matrix whose diagonal and upper right elements 
are independently and identieally distributed random variables with mean 0 and varianee 
1. Then it follows that eonverges in the sense of probability to psc- Let R denote 

the variable spaee for the i.i.d. random variables, then 


lim inf P 

n—>-oo 


(l£ 


hiy)dpB„/v^iy)- h{y)dij,sc{y) 


/ X 



1 


for all e > 0 and every bounded eontinuous funetion h : R ^ R. 

A conservative version of [5, Theorem A] is now stated 

Theorem 3. Let Mn he a symmetrie nxn matrix whose upper right elements are inde¬ 
pendently and identieally distributed random variables with mean 0 and varianee , and 
whose diagonal elements are independently and identieally distributed with mean 0 and 
finite varianee. Then it follows that \Xi(Mn)\ = 2cry^(l + o(l)) asymptotieally 

almost surely. Stated another way. 


lim inf P ( sup |Ai (Mn) | = 2cr\/n(l + o(l)) )= 1. 
n^oo \i<i<n J 

Theorem 1 states that the spectral distribution of a random matrix with elements 
drawn from a normal distribution with mean 0 and variance 1/n converges to the unit 
disk centered at the origin of the complex plain. For the same random matrix but with 
the added assumption of symmetry, the spectrum converges to the line segment [—2,2] 
on the real line, see Figure T3. The spreading of the spectrum from a diameter of 2 to 
a diameter of 4 can be explained by the fact that all 2-cycles in the matrix now have a 
positive loop game. Positive feedback for 2-cycles always repels the eigenvalues along the 
real axis [30]. 


4 Stability 

4.1 Primer on Stability 

Consider the time-varying dynamical system defined by 

x{t) = f{x{t),t) 
x(to) = Xq 

where x : R ^ R’^ is the state vector, t is time, to is the initial time, and ( ) — ^( )• 
Let X* G R’^ be the equilibrium solution so that /(x*,t) = 0 for all t. The solution to 
the ordinary differential equation in (T7) is a transition function 0(t;xo,to) such that 
(/)(to; xo, to) = xo and 0(t;xo,to) = /(0(t; xo, to), t). For existence and uniqueness condi¬ 
tions see [16]. Below we give the various definitions of stability as defined in [40,48,63,70]. 
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[A]ij ~ Af{0, 1)/[A]y = [A]ji ~ 7V(0, l)/v^ 



Real Real 


Supplementary Text Figure T3: Eigenvalues for (left) random matrix and 
(right) symmetric random matrix of dimension n X n with n = 1000 where ele¬ 
ments are drawn from the distribution A/'(0, V)l\/n. 


Definition 2 (Stability). Let to > 0, the equilibrium x* is 

(i) Stable, if for all e > 0 there exists a S{e,to) > 0 such that ||a:o — x*\\ < S implies 

|| 0 (to; to, ^co) — < e for all t > to- 

(ii) Attracting, if there exists a p(to) > 0 such that for all 77 > 0 there exists a T(irf, xq, to) 

such that ||xo — || < p implies || 0 (t; xo, to) — ai* || <77 for all t > to + T. 

(iii) Uniformly Stable, if the S in (i) is uniform in to, thus taking the form (5(e). 

(iv) Uniformly Attracting, if it is attracting where p does not depend on to and T(rj,p) 
does not depend on xq or to- 

(v) Uniformly Asymptotically Stable (UAS), if it is uniformly stable and uniformly at¬ 
tracting. 

(vi) Uniformly Bounded, if for all r > 0 there exists a B(r) such that ||xo — ^ t 

implies that ||s(t;to,xo) — < B for all t > to- 

(vii) Uniformly Attracting in the Large, if for all p > 0 and 77 > 0 there exists a T(p,p) 

such that ||xo — || < p implies ||s(t; xq, to) — a;* || <77 for all t > to + T. 

(viii) Uniformly Asymptotically Stable in the Large (UASL), if it is uniformly stable, uni¬ 
formly bounded, and uniformly attracting in the large. 

(ix) UAS in the Positive Orthant, if it is uniformly stable, uniformly bounded, and uni¬ 
formly attracting in the positive orthant. 

The precise definition of UAS is important in the context of dynamical systems. If one 
is able to show that a given system is UAS, then it follows that the dynamics are stable in 
the presence of bounded disturbances as well. That is, if the dynamics in (T7) are UAS 
then for ||d(t)|| < 7 where 7 > 0 is sufficiently small, the dynamics 

+ d{t) 

are uniformly stable. If the dynamics in (T7) are UASL then x{t) is bounded for all 
bounded d{t) and 7 can be arbitrarily large. A detailed discussion regarding this fact can 
be found in [40, Definition 56.1 and Theorem 56.4] and a practical example in the context 
of adaptive systems can be found in [69]. 

For a linear dynamical system xft) = Ax(t) the uniform asymptotic stability (which is 
actually exponential) is verified if all of the eigenvalues of A have real parts less than zero. 
A well known theorem regarding the stability of linear systems, due to Lyapunov, is now 
stated. 

Theorem 4 (Lyapunov). The eigenvalues of a real matrix A have all real parts less than 
zero if and only if there exists a P = > 0 such that A^P + PA < 0. 
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A stronger version of Thereom 4 will be needed when discussing the stability of Lotka- 
Volterra dynamics. 

Definition 3. If there exists a diagonal positive matrix P such that A^P + PA < 0 then 
A is said to be Diagonally Stable. 

4.2 Stability of Generalized Lotka Volterra Dynamics 

Consider dynamics of the form 

n 

Xi{t) = riXi{t) P Xi{t)'^aijXj(t), z = (T8) 

i=i 

where t E [to,oo) is time with to the initial time. The state vector is denoted x E 
and defined as x = [xi, X 2 , ..., The linear terms are collected in the column vector 

r = [ri, X 2 , Vn]^ and A = (cLij)i>i,j>n captures the pair-wise interactions in the 
generalized Lotka-Volterra dynamics presented in (T8). The dynamics in (T8) can be 
compactly represented as 

x{t) = diag(x(t))(r + Ax{t)). (T9) 

A discussion regarding the invert ability of A is in order. This will become important in 
determining whether the system in (T8) has a unique non-trivial steady state. As discussed 
in the main text, the existence of the Verhulst terms increases the likelihood that A 
is full rank. First consider the case where all aij = 0 for all i ^ j, then a necessary and 
sufficient condition for A to be full rank is that all diagonal elements are non-zero. 

If A is invertible then there exists a unique non-trivial steady state solution of the 
dynamics in (T8), denoted as x* = —A~^r [38]. We are interested in answering the follow¬ 
ing question. Given an ecological system of n species, is it possible to introduce another 
species and drive the system to any non-trivial steady state of our choosing? We will show 
that the answer is yes, and furthermore, if the original n-dimensional ecological system 
satishes a diagonal stability condition, then we can design the interaction strengths for the 
introduced species so that the new n + 1 dimensional ecological system is asymptotically 
stable for all initial conditions in the positive orthant. We then discuss equivalent results 
regarding steady state shift and stability when an arbitrary number of species is added. 
Consider the following set of assumptions. 

Assumption 1. A is invertible. 

Assumption 2. For the dynamics in (T8) the steady state solution x* G M^o- 
Assumption 3. The matrix A is diagonally stable, see Definition 3. 

Theorem 5 ([38, Theorem 1]). If the system in (T8) satisfies Assumptions 1-3, then the 
steady state x* is uniformly asymptotieally stable for all initial eonditions xo G M>o- 

Proof. Let V (x, t) = 2 Pi{xi — x* — x* log(xi/x*)) be the Lyapunov candidate where 
Pi is the z-th diagonal element of a diagonal positive matrix P such that A^P + PA < 0. 
Differentiating the Lyapunov candidate it follows that 

V{x,t) = —) 

1 \ / 

1 = 1 ^ ^ 

= 2'Upt{xt - X*) — 

U 

n / n 

= 2^Pi(xi - X*) I bi p'^aijXj 
i=l \ j=i 

n n 

= ‘2y2pi{xi- X*) aij(xj - X*) 
i=l 3=1 

= (x- x*y(A'^P + FA)(x - X*). 


Thus the Lyapunov candidate is positive dehnite in x — x* and its derivative is negative 
definite in x — x*. □ 
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Remark 1. Note that in the original work of Goh [38], stability in the positive orthant 
is denoted, but in fact what he proved is uniform asymptotic stability in the positive 
orthant. It is important to distinguish the two as stability only implies boundedness of 
trajectories and uniform asymptotic stability implies convergence to the equilibrium as 
well as robustness to bounded persistent disturbances. 

4.3 Stability in the presence of new species 

4.3.1 Adding one new species 

For the following examples we are interested in the m = n + 1 dimensional dynamics 

z(t) = diag(z(t))(s + Fz{t)) (TIO) 


where 


r 

, and F = 

A b 



s 


c d 


with A and r are as dehned in (T8), and we have introduced the new elements s,d G M, 
and b,c eMT. These dynamics represent the addition of one new species to the ecological 
system in (T8). 

We now introduce an important property of diagonally stable matrices. 


Theorem 6 ([84, Theorem 3.1]). Let A G be partitioned as 


All Ai2 
A21 A22 


where An G ^ ^12,^21 G and A 22 < 0. Then A is diagonally stable if 

and only if An and the quantity (An — A 12 A 21 /A 22 ) have a common diagonal Lyapunov 
function, i.e., there exists a positive diagonal P such that AjiP + PAn < 0 and {An — 
Ai 2 A 2 ilA 22 yP + P{An — A 12 A 21 IA 22 ) < 0. 

Lemma 1. For any steady state solution x* of (T8) satisfying —Ax* = r there exists 
b,c,d,s such that any z* G M>o can be made to be a steady state solution of (TIO). 


Proof. Any steady solution of (TIO) satisfies the relation g = —Fz*, which when expanded 
denotes the following relation 


r 

s 


'A b 


A* 

^l:n 

d 


0 


There are 2n + 2 degrees of freedom in the variables b, c, d, s and there are only n + 1 
constraints. The variable b is fixed by the top row of the above equation and can be 
expressed in closed form as 

(Til) 

Then for any c and d, s can be chosen as 


S = —c Zi-. 


— dz^ 


□ (T12) 


Corollary 7. For any steady state solution x* of (T8) satisfying —Ax* = r, and given 
any c G 'MT and d G M there exists b, s such that any z* G M>o can be made to be a steady 
state solution of (TIO). 

Theorem 8. For the dynamics in (T8) with p = 1 satisfying Assumption 3 there exists 
b, c, d, s such that any z* G M>o can be made to be asymptotically stable for all initial 
conditions in the positive orthant. 


Proof. From Lemma 1 it follows that for any z*,c and d the elements s,b are fixed. Now 
we will show that for any c there exists a d such that z* is asymptotically stable. We begin 
by assuming that b and s are fixed by (Til) and (T12). From Assumption 3 we know that 
A is diagonally stable, and thus there exists an e > 0 and diagonal matrix P > 0 such that 
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A^P + PA < — e/, where I is an appropriately dimensioned identity matrix. Choosing 
d < 0 and |d| > 2Xma.^{P)\\bc^\\/e, and noting that 

(A - bc^/dyP T PiA - bc^/d) < -elT2X^,yP)\\bc^\\/d 

we can deduce that {A — bc^/d^P + P{A — bc^/d) < 0. Thus {A — bc^jd) and A have 
a common diagonal Lyapunov function. This, in addition with the fact that d < 0, 
the conditions of Theorem 6 are satisfied and thus F is diagonally stable. Now, applying 
Theorem 5 to the dynamics in (TIO) we deduce that z{t) can be made to be asymptotically 
stable for all z(to) G □ 

4.3.2 Adding an arbitrary number of species 

For the following examples we are interested in the m = n + p dimensional dynamics 

z{t) = dmg{z{t)){g + Fz{t)) (T13) 


where 



and 



B 

D 


with A and r are as defined in (T8), and we have introduced the new elements s G 
D G and B,C E These dynamics represent the addition of p new species to 

the ecological system in (T8). 


Theorem 9. For any steady state solution x* of (T8) satisfying —Ax* = r there exists 
B,C,D,s sueh that any z* G R>o can be made to be a steady state solution of (T13). 
Furthermore, if A is diagonally stable, then B,C,D,s ean be ehosen sueh that the system 
in (T13) is uniformly asymptotieally stable in the positive orthant. 


Proof. Any steady solution of of (T13) satisfies the relation g = —Fz* which when ex¬ 
panded denotes the following relation 


r 


A 

B' 


^l:n 

s 



D 




There are 2np + + p degrees of freedom in the variables B,C,D,s. The variable B is 

fixed by the top row of the above equation and must satisfy the following relation 




There are n x p degrees of freedom in the selection of B and only n constraints in the 
above equation. Thus such a B always exists for any r. A, and z*. For any C and D, s 
can be chosen as 

S C Zi-^ 

Finally, we show that with the extra degrees of freedom in C and D there always exists a 
diagonal Pi > 0 and P 2 > 0 such 


'A 

B' 

T 

Pi 

o' 


■Pi 

o' 

'A 

B' 


D 


0 

^ 2 . 

+ 

_0 

^ 2 . 


D 


(T14) 


Given that by assumption there exists a diagonal Pi > 0 such A = A^Pi + Pi A, then by 
the Schur complement the inequality in (T14) holds if an only if 

D^P 2 + 2 P 2 D - + P 2 C^)A~\CP 2 + PiB)< 0. 


Given any A, B,C and positive diagonal Pi , P 2 there always exists a D such that the above 
inequality holds. □ 
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4.3.3 Removing an arbitrary number of species. 

Diagonally stable matrices have a very special property that every principle minor is also 
diagonally stable, giving us the following definition and theorem. 

Definition 4. Let L be a proper subset of N = {1, 2,... ,n} then a principle minor of 
A G is the matrix obtained by omitting the columns and rows of A whose index 

appear in L. 

Theorem 10 ([19, Theorem 1]). If A ^ R^^’^ is diagonally stable, then all principle 
minors of A are also diagonally stable. 

In the context of Lotka-Volterra dynamics this implies that if a given systems is diagonally 
stable, then even if an arbitrary numbers of species are removed from the system, the 
resulting system is still uniformly asymptotically stable in the positive orthant. 

Corollary 11. A necessary condition that A is diagonally stable is that all of its diagonal 
elements must be strictly negative. 

4.4 Robustness to disturbances. 

We now consider three classes of Lotka-Volterra dynamics in the presence of disturbances, 
the first two are determinstic and the second is stochastic. 

4.4.1 Deterministic Dynamics 

In this section we will analyze the following two dynamical systems 

x{t) — diag(x(t))(r + Ax{t) + w(t)) (T15) 

and 

x(t) — d(t) + diag(x(t))(r + Ax{t)) (T16) 

where d{t) and w{t) are known to be a priori bounded. In (T15) the term w{t) represents 
uncertainty in the growth rates of the species. While in (T15) the term w is deterministic, 
in the since that it is bounded, in terms of the ecology literature this term is sometimes 
referred to as the stochastic effect. Note that we will treat this term in a purely stochastic 
since shortly, complete with Ito calculus. The term d{t) is a migration term. The following 
theorems show that in the presence of these disturbances, the systems are bounded for all 
initial conditions (in the positive orthant). 

Theorem 12. If the system in (TI5) satisfies Assumptions 1-3 with ||d(t)|| < a then the 
state X is uniformly bounded for all initial conditions xq G R^o- 

Proof. The proof is given for two different scenarios. The first scenario (when the distur¬ 
bance is small) uses the same Lyapunov function as was used in the analysis of the distur¬ 
bance free dynamics. In the second scenario a new Lyapunov candidate is introduced. A 
few definitions are needed before the different scenarios are analyzed. Let P be a diagonal 
positive solution to the Lyapunov equation A^P + PA = — Q, where Q = > 0, and 

Pi denotes the z-th diagonal element of P, just as in the disturbance free case. Let ^min 
denote the minimum eigenvalue of Q and Pmax denote the maximum diagonal element in 
P. Scenario: (a) the compact set 

{x : \\x-x*\\ < 2p max Q^/^min} 

does not intersect any of the n-axes; and (b) when the above compact set does intersect 
at least one of the n-axes. 

We now address the stability of Scenario (a). Let 

n 

V{x) = 2'^pi{xi - X* -X* log(a:i/a:*)) 

i=l 


(T17) 
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be a Lyapunov candidate. Differentiating the Lyapunov candidate it follows that 

= ‘^Ypi (xi - XiE\ 

1 \ / 

1=1 ^ ^ 

= 2'^pi{xi - X*) — 

. Xi 

Z =1 

n / n 

= 2'^Pi{xi - Xi) I bi + di{t) + y^^ajjXj 

i=l \ 3=1 

n / n 

= 2^^Pi{xi - X*) I di(t) + Y^aij(xj - X*) 
i=l \ 3=1 

= (x — x*)^(A^P + PA){x — a:*) + 2{x — 

Recall that ^min denotes the minimum eigenvalue of Q and Pmax denotes the maximum 
diagonal element in P, then it follows that 

V < -qmirv\\x - + 2pmax||iC “ X^\\a. 

Thus for all ||a; —x*|| > 2pmaxQ^/^min, R < 0 and thus V{x(t)) < oo for all t > to- It follows 
that for all xo G M>o the state x(t) is bounded. Note that x{t) can never be negative due to 
the fact that the disturbance appears as Xi{t)di{t). This completes the proof for Scenario 
(a). 

Scenario (b) has to be treated differently do to the fact that if any of the = 0, then 
it follows that for the Lyapunov candidate in (T17) V = oo. This was not possible in 
Scenario (a), but is possible in Scenario (b). We dehne a new Lyapunov candidate 


y(x) = 2y^Vi(xi) 

i=l 


(T18) 


where 


Vi{x) 


Pi{Xi - X* - X* log(Xi/x*)) X* < Xi 

0 0 < Xi < X*. 


For a general set of dynamics the above candidate function could not be a Lyapunov 
candidate. However, do to the special form of the Lotka-Volterra dynamics x(t) G P>o for 
all time regardless of the specihc coefficients. Any population dynamic model should have 
this property, as negative abundances would make no sense. Differentiating V in (T18) we 
have that 

n 

V{x) = 2 y^Ui(xi) 

i=l 


where 


Vi(x) = 


Pi{xi - Xi) (di{t) + aij{xj - X*)^ 

0 


Xi < Xi 

0 < Xi < Xi, 


which is continuous in x, note that pi{xi — x*) \^di{t) T ciij{xj — x*)j = 0 when 

Xi = x*. 

We now dehne the index set for all species as 


X = {z : 0 < z < n, z G N} (T19) 

and the set of indices for species abundances that are greater than x* is dehned as 

<S = {z : Xi > X*} C X. 

Using this notation we can write the derivative of the Lyapunov function as 


V{x) 


2'^Pi{xi - x*) 
ies 


( di + aij{xj -x*)+J 2 I 

V jes jes^ / 


(T20) 
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The above equation can be rewritten in a more compact matrix-vector form as 

V(x) = (xs — xs)^ (As,sJ^s,s + Ps,sAs,s ] (xs — x^) + 2(xs — xs)^Ps,sds 

^ ^ (T21) 

+ 2 (X5 — Ps,sAs,S^ (xs^ — Xsc) . 

Using Theorem 10 we can exploit a very unique property of diagonal Lyapunov functions, 
and that is the diagonal stability of all principle submatrices. It is also easy to confirm that 
the same P matrix can be used for the sub-diagonal Lyapunov equations. Therefore, it 
follows that A^ ^Ps,s T Ps,sAs,s = —Qs,s < 0 for any S. Using the following definitions 

= minAi(Q5,5) 

I 

Pmax(>5) = maxAi(Ps,s) (T22) 

I 

p\S) = \\AsM\, 

the equality in (T21) can be bounded as 

V{x) < -gmin(‘5)l|a:5 + 2pmax(<5)||a;5 - XsW (a+ ^'(5) 11*5' -®s^ll) • 

Given our definition of S it follows that \\xs^ — < \\x%c\\ < ||x*||. This allows the 

above inequality to be further simplified as 

Vix) < -gmin(>5)||a:5 - xsW^ + 2pmax(5)||x5 “ 2:51| (a + /?'(5)||a:*||) . (T23) 

The bound in (T23) is still not sufficient to address stability as the bounds depend on the 
set of indices in <S. In a two step process we obtain uniform constants and then introduce 
a distance function that is independent of S. We will now give uniform analogs to the 
bounds ^min(‘^)?Pmax(<5)?by taking the appropriate supremum or infimum over the 
powerset 2^ of all subsets S C X. The powerset of interest is finite and thus the following 
are well defined 


q = 


inf 

5CX\0 


q'miniS) 


p= Suppmax(>5) 

sex 


p = sup ^'(5). 

sex 


(T24) 


Using the uniform bounds above the inequality in (T23) can be replaced by 

y{x) < -q\\xs - xsW^ + 2p\\xs - xJII (a + ^||a:*||) . (T25) 

We now define a distance metric between an element y E MX and a compact set ^ G 
as 

d{y,A) = \ni{\\y - z\\ ■. z e A}. (T26) 

Given that A is compact, such o, z ^ A always exists. Our compact set of interest is 
defined as 

X = {x:^<Xi< X-}. (T27) 

Using the compact set defined just above, the inequality in (T25) can be equivalently 
rewritten as 

V (x) < —qd(x^ X)^ + 2pd{x, X) (a + ||) . 

Thus for all d{x,X) > 2p (a + ^||x* ||) /q, U < 0. Therefore V{x{t)) is bounded for all 
t > to. □ 


Theorem 13. If the system in (T16) satisfies Assumptions 1-3 with ||d(t)|| < a and 
furthermore we exelude the possibility for the disturhanee to generate negative state values^, 
then the state x is uniformly hounded for all initial eonditions xo G M>o- 

^This is not a restrictive assumption in the context of ecology as a disturbance can never result in 
the creation of a negative population. 
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Proof. Consider the Lyapunov candidate in (T18) and differentiating along the system 
dynamics in (T16) it follows that 

V{x) = 2'^pi{xi - Xi) ( di/xi + '^aij(xj - Xj) + ^ CLij{xj - Xj) \ . 
ies \ jes jes^ ) 

where S C I was defined in (T19) and (T20). Rearranging terms slightly with regard to 
di /Xi it follows that 

=2'^pi{xi - X*) ('^aij(xj -x*)+ aij(xj - x*) | + 2'^diPi{l -x*/xi). 
ies \jes jes^ / ies 

Given our definition of S it follows that (1 — x'l/Xi) < 1 when i E S. Therefore it follows 
that 


V(x) < 2y^pi(xi 
ies 


jes jes^ / ies 


The above inequality can be rewritten in a more compact matrix-vector form as 

V{x) < {xs — (^^s,sPs,s + Ps,sAs,s^ (xs — xs) + 2||P5,51| || 

+ 2(xs - xs)^ Ps,sAs,s^ (xs^ - xsc). 

Using the bounds in (T22) the above inequality can be reduced to 

V(x) < -qLin(S)llxs - + 2p'^^^(S)llxs - x*s\\f3'{S)\\xs- - xJcH + 2p^^^{S)a. 

Given our definition of S it follows that Ijx^c — x%c\\ < H^JcH < \\x*\\. This allows the 
above inequality to be further simplified as 

V{x) < -^min(<S)||x5 - + 2p^ax(<S)||x5 “ 11(<S)11X* 11 + 2p^^^{S)a. (T28) 

Using the definitions in (T24) and the set distance function in (T26) with the compact set 
X defined in (T27) the above inequality can be written as 

V (a:) < —qd{x, X)^ + 2pd{x, X)^\\x* || + 2pa 

Rewriting the above inequality as 

V{x) < -|d(x, X)‘^ - I ^d{x, X)‘^ - 4^d{x, X)0\\x*\\^ + 2pa 

and completing the square with respect to the middle term | (^d(x, X)^ — 4^d(x, X)^\\x* ||^ 
the following inequality holds 


V{x)<-Y{x,Xf-^\^d{x,X)-2^^x*\\^ +2pa + 2 

Noting that the term —| {d{x^ X) — 2|^||x*||^ < 0 it follows that 

V{x) < -U{x,Xf+ 2 pa + 2 ?AlEL 
2 q 

From the inequality in (T29) it follows that 


Paw 


2 /Q 2 II ||2 


,rn\x 


(T29) 


d{x,X)>2J^ + P"^Y^^" 

V 9 

implies U < 0. Therefore, x{t) is uniformly bounded and asymptotically converges to the 
compact set 




pa 
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4.4.2 Stochastic Dynamics 

In this section we wish to analyze a stochastic differential equation that is similar to (T15), 
but we no longer make the assumption that the disturbance is bounded, instead we assume 
that the disturbance appears as a Brownian motion w(t) G MT, resulting in the following 
differential equation, 


dx = diag(a:)(r + Ax) dt + c diag(x) dw. (T30) 

The variable c G M will be a constant used to scale the square root of the variance of the 
brownian motion. It will be shown that x(t) converges in a probabilistic since to a compact 
set which is proportional to c. The analysis of stochastic differentials is signihcantly more 
challenging than their deterministic counterparts. Before stating the stochastic version 
of the stability result, we first need to formally dehne filtration. Brownian motion, and 
give a key result of Ito. In this work the linear ordered set of interest will always be time 
t G [to, oc). All of the following definitions are given in a less general context compared to 
their original presentation [24]. 

Definition 5 (Stochastic Process, [24, Part 2 §1.8]). Let (Q, P) be a probability space. 
A stochastic process on (Q, T, P, [to, cxo)) is a family of random variables {y(t), t G [to, oo)} 
with map (t,ct;) yit^u) from [to,oo) x Q to it, the variable space of interest. Given 
that each instance in time y{t) is a random variable, the variable space R = {R,R) is 
necessarily a measurable space. 

Definition 6 (Filtration, [24, Part 2 §1.1]). Let (12, R) be a measurable space and t G 
[to, oo) is time, then a filtration of that measurable space is a map t J^(t) of increasing 
sub a algebras such that R(s) C R(t) C R when s < t. 

Definition 7 (Adaptation, [24, Part 2 §1.1]). Let {y{t),t G [to, oo)} be a stochastic process 
from a filtered measure space {Q, J^(')} into the measurable space {R,IZ}. The process 

is adapted to R(‘) if for each t, y(t) : (LL,R(t)) (R,R) is measurable. 

Definition 8 (Progressively Measurable, [24, Part 2 §1.2]). Let {y{t), R{t),t G [to,oo)} 
be an adapted stochastic process. The function y : [to,oo) x 12 ^ R, where (R,1Z) is a 
measurable space, is deemed progressively measurable if 

y : ([to, t] X Q, Borel([to, t]) x R(t)) (R, IZ) 

is measurable for any t > to. 

Let us pause now and discuss these dehnitions. Dehning a stochastic processes is rather 
obvious. We extended the definition of a random variable to incorporate time. At each 
hxed instance a stochastic process is nothing but a random variable. The extra dehnitions 
and the progression from hltrations, adaptations, and progressive measurability, are needed 
so that we can better understand what is needed to perform the following integration 
[88, Remark 7.1.1]. Let {y(-)?‘^(')} be an adapted stochastic process dehned as before 
where y : [to,oo) x 12 ^ R. Let f : R ^ R he bounded and IZ measurable, then the 
map (t,uj) z — f(y(r,(ju)) dr need not be adapted. However if in addition to being 

adapted the stochastic process is progressively measurable, then z will be progressively 
measurable as well. This process is a key necessary ingredient when studying stochastic 
differential equations as it can be inherited. Next we formally dehne a Brownian motion 
and give a sufficient condition for the adaptation so that the progressive measurability 
condition is satished. 

Definition 9 (Brownian, [24, Part 2 §VIL2 ]). A Brownian motion is an adapted stochastic 
process {le(-), J^(-)} from the hltered probability space (Lt,R,B^R(t),t G [to,oo)) into the 
measurable space (RT, Borel(R’^)) 

• that is Markovian, i.e. when s < t and A G Borel(R’^) it follows that 
B(w{t) G M|J’(s)) = P(^(t) G A\w{s)) 


almost surely. 
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• with a stationary stochastic transition function (the transition function is indepen¬ 
dent of the current time of the process) with density — ^) defined on R x R’^ 

relative to the n-dimensional Lebesgue measure, where p is dehned as 


p{t,'ip;a, to) 


(27rcr^t) ^^^exp if t > to 

0 if t < to, 


where t G R is time and '0 G R’^ is the space variable. 


• that is almost surely continuous 

We are almost ready to address the problem at hand, however we still need to show 
how we can ensure that the Brownian motion is progressively measurable. Unfortunately 
we need a few more definitions and then we will show that Brownian motions can be 
naturally adapted so as to have the progressive measurability property. 

Definition 10 (Right Continuous Filtration, [24, Part 2 §1.1 ]). Let (Q, T, T(t)U ^ 
[to,oc)) be a filtered measurable space, and define — C\s^tT{s) for all t G [to,oo). 

The filtration is right-continuous if J^(-) = *F^('). 

Theorem 14. Let {y(-)?*^(‘)} ^ Brownian motion into R’^ with respect to time t G 

[to,oo), and if T{t) is generated by the null sets and a{y{s);s < t), the smallest sigma 
algebra for which all y{s) are measureable for all s G [to, t], then J^(-) = -^^(O 

Proof, see [24, Part 2, §VI, Theorem 8] □ 

This theorem implies that one can assume without loss of generality that the filtration 
is right-continuous. That is, there is a natural way to construct them for any brownian 
motion. 

Theorem 15. Let {y(-)?*^(‘)} ^ Brownian motion with right-continuous filtration con¬ 

taining the null sets, then the Brownian motion is progressively measurable. 


Proof, see [66, Part A Theorem 47] 


□ 


As stated before, this progressive measurability is necessary when discussing the so¬ 
lutions to stochastic differential equations as it allows for integrals containing stochastic 
processes to inherit the progressively measurable property. We now state a classic result 
do to Ito [42]. Our version follows from [24, Part 2 §VIII.12]. 

Lemma 2. Let {'^^(•)?^ Brownian motion into R’^ with a right-continuous filtra¬ 
tion containing the null sets. Consider the dynamics x{t) G R’^ given by 

dx = pdt + adw 

where t G [to, oo), p{x) G R’^ is locally Lipschitz in x, and a(x) G R^^’^ is globally Lipschitz 
in X. Assuming (x, t) i— f{x,t) G R ts twice differentiable and continuous in terms of x 
(class C^ with respect to x) and once differentiable and continuous in terms oft (class C^ 
with respect to t), then 

d/= + (V./)Gdt + ^ tr((7Hess.(/)(7) dt + (V./)Vdw. 

Remark 2. In most constructions it is assumed that p is globally Lipschitz. Indeed this 
assumption was made in one of Ito’s original papers [43]. A detailed discussion regarding 
the existence of solutions when only locally Lipschitz conditions are assumed can be found 
in [78]. Exploiting the stability of our dynamics the existence and uniqueness of solutions 
can be deduced directly from the original work of Ito however. In our analysis we will show 
that for the dynamics of interest the state is uniformly bounded almost surely, a formal 
definition is given just below this remark. Therefore, we can define equivalent dynamics 

dx — pdt + adw 

where p = p on the compact set of interest, and zero outside this compact set. Given that 
p is locally Lipschitz it follows that the function p with compact support is by dehnition 
globally Lipschitz. Solutions x exist, are unique, and are uniformly bounded almost surely. 
This then implies the existence and uniqueness of x for the dynamics dx = pdt + adw on 
the same compact set which x(-) remains in, almost surely. 
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Definition 11. The state x(t;to,xo) as a solution to the difference equation in (T30) is 
uniformly hounded with probability one if for ever r > 0 there exists a B(r) >0 such that 
||^(^o)|| < r implies 

P ^sup||a:(i)|| < B(r)j = 1 

for all t > to. An equivalent statements would be x(t) is uniformly bounded almost surely. 

Theorem 16. Consider the dynamies in (T30) for x{t) G MT. Let assumptions 1-3 hold 
with {re(-), a Brownian motion into BT with a right-eontinuous filtration eontaining 

the null set. The following then hold 

1. The state variable x(-) is uniformly bounded in expeetation, 

2. x{-) is uniformly hounded with probability one, and 

3. x[‘) asymptotieally eonverges to the eompaet set 

V ^ [x : \\x - x*\\ < c.TTLTLTX (t 31 ) 

[ V ^min J 

with probability one. Stated more preeisely, 

P { lim x{t) G p) = 1. (T32) 


Sketeh of the proof. This proof only outlines the analysis when the compact set V does 
not intersect any of the n-axes. The more complicated scenario when this is not the 
case can be handled just as it was in the proof of Theorem 12. Consider the Lyapunov 
candidate V{x) = — x* — x* \og{xi/xi)). Taking the differential along the 

lines of Lemma 2 (Ito) results in 


dV = {\/xV)^dmg{x){r + Ax) dt + ^ tr(diag(x)HesSx(y) diag(x)) dt 

+ (Va.y)’^c diag(x) dw (T33) 

Recall from the steps in the proof of Theorem 12 that 

(VccR)^ diag(a:)(r + Ax) = —{x — x*)^ Q{x — x*) 
where A^P + PA — —Q. Thus giving 


dV = —(x — — X*) dt + ^ tr(diag(x)HesSa;(y) diag(x)) dt 

+ (Va.R)’^c diag(x) dw (T34) 


Next, note that [HesSa;(R)]zz = 2^^ and [HesSa;(R)]zj = 0 when i ^ j. Substitution into 
the above equation results in 

dV — —{x — x*)^Q{x — X*) dt + dt + (VxV)^c diag(x) die (T35) 

i 

Note that x{t) and die(t) are independent, x{t) is only dependent on die(s) when s < t. 
Therefore E((Vccy)^c diag(x) die) = 0, given that for a Brownian motion Edie = 0.^ 
Therefore 

< -^min||x - +C^nXmaxPmax (T36) 

where ^min is the minimum eigenvalue of Q, Pmax is the maximum diagonal element in P, 
and Xmax is the value in the vector x*. Therefore, 


\\x — X W > c 


^maxl^max 


^min 


EdR 

dt 


< 0. 


(T37) 


^In order to address this rigorously we would rewrite the expression in (T35) in integral form. Then 
we would integrate only up to a stopping time which we designed so that {pic. diag(a::) < oo up to 
the time of interest. Then we would use the fact that Yigw = 0 for any bounded g. The bound used to 
generate the stopping time would then be relaxed and the desired result would be obtained as t ^ oo 
using Fatou’s Lemma and monotone convergence. A similar procedure is carried out in [21, (4.3)-(4.5)] 
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Given that EdV/dt < 0 outside a compact set it follows that EV is uniformly bounded in 
terms of the initial condition x(to) [50, Lemma 5.4]. Given that V{x) > 0 is continuous 
in X and convex it follows from Jensen’s inequality that Ex is uniformly bounded as well 
[47]. Jensen’s inequality will be used without reference from this point forward. Glaim 1 
has been proven. 

We now move on to the second claim, namely that the dynamics are uniformly bounded 
with probability one. We address this claim with sub scenarios: (a) x{to) ^ P, and (b) 
x{to) E V. Under scenario (a) we will now show that P(supU(x(t)) — V{x{tQ) < Ji) = 1 
for any Ji G (0,oo). The previous statement is equivalent to P(supU(x(t)) — U(x(to) > 
Ji) = 0. This statement will be proved by contradiction. Assume that there exists a 
^2 E (0,1] such that P{supV{x{t)) — V{x{tQ) > Si) > 82 - Using Markov’s inequality it 
follows that 

EsupU(x(t)) — U(x(to)) > 8182 - 

This statement however contradicts (T37) which states that outside the compact set of 
interest the expected value of the Lyapunov candidate is strictly decreasing. Therefore 
it follows that P{supV{x{t)) — V{x{to) > Ji) = 0 which is equivalent to the claim that 
P(supU(a:(t)) — V{x{to) < Ji) = 1. Application of Jensen’s inequality then implies that 
X is uniformly bounded with probability one for sub-scenario (a). 

We now address sub-scenario (b) where it is assumed that x{to) E V. If a:(-) remains 
in P for all time, then we are done and the state is uniformly bounded. Therefore, assume 
at some time ti the state x{ti) ^ V. This sub-scenario is now equivalent to sub-scenario 
(a) with time shifted. Therefore, it follows that x is uniformly bounded with probability 
1 in sub-scenario (b) as well. This completes the proof of claim 2. 

We approach claim 3 following a method proposed in the proof of [21, Theorem 2.1]. 
Gonsider the probability of three mutually exclusive scenarios, just as in [21, Theorem 2.1] 


Pi = P [ limsupd(x(t), P) = 0 | 

\ t^oo J 

P 2 — P fliminf d(x(t),P) > o) 

V t^oo J 

P 3 = P (liminf d(x(t), P) = 0 and limsupd(x(t), P) > 0 ) 

V t^oo J 

We wish to prove that pi = 1 and P 2 ,P 3 = 0. 

We hrst prove that p 2 — 0. This will be proved by contradiction. Assume p 2 = ei 
where ei G (0,1], then there exists an €2 such that P (liminU^oo d(x(t), P) > 62 ) = ei. 
Using Markov’s inequality it follows that 

— E lim inf d(x(t), P) > ei. 

62 t^oo 

Multiplying both sides by 62 it follows that 


Eliminf d(x(t),P) > 6162 ■ 

t^oo 

From the dehnition of liminf it follows that for any 63 G ( 0 , 6162 ) there exist a hnite 
Ti G [to, 00 ) such that 

Ed{x{t),V) > 63 for all t > Ti. (T38) 

Therefore, for t > Ti it follows that EdV/dt < 0. Given that V{x*) = 0 and V{x) is strictly 
increasing away from x — x* and tends to inhnity for large x in the positive orthant it 
follows that E||x — X* II is strictly decreasing in time. Given that P is centered at x*, there 
exists T 2 such that 

Ed{x{t)^P) < 63 for all t > T 2 . 
which contradicts (T38). Therefore p 2 — 0. That is 

P (liin ini d{x{t)^P) > o) = 0. 

V t^oo J 

We now establish that p 3 = 0. As before we achieve this by contradiction. For any 
65 > 0 assume 


(liminf d(x(t),P) = 0 and limsupd(x(t),P) > 65 ) 7 ^ 0 . 


V ^ 


(T39) 
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For any ee and er such that 0 < er < ee < € 5 , consider the stopping times'^ T[ \ LL ^ [to, 00 ), 
T/' : Q ^ [to, 00 ) with i G N>o, where 

T' = inf {t > T/ii : d{x{t),V) < 67 } 

T'' = inf {t>T' : d{x{t),V) > ee} . 

From (T39) and the stopping times dehned above, it follows that 

lim T/ = 00 and lim T-' = oc 

z—>-oo z—>-oo 

In order to continue with this proof by contradiction we need to obtain a lower bound on 
the expectation T/^i — T/^ Indeed if this can be done, then we expect the solutions to the 
differential equation to spend a non negligible amount of time within a domain of the state 
space where the Lyapunov function will be decreasing. Thus we can see how this may lead 
to a contradiction, for how can trajectories be expected to spend an inhnite amount of 
time in a domain where the Lyapunov function is always decreasing. Following the same 
procedures as in the proof of the bound in [21, (2.27)] it can be shown that there exists 
an 68 > 0 such that 

B{T'^,-T'' \T{T''))>e8. 

Using the bound above and the definition of V in (T31) it follows that 


E ^y* qmin\\x{T) - X^'lf 


dr T'{T'') 


> 


^^maxPmax 


I ^ 

^6 ~r ^minC6 


^8 


From the bound in (T36) it follows that 

POO 

V(x{to)) >E / qniin\\x{T) - - C^na^maxPmax dr 

JtQ 

( I ^^maxPmax 2^ Tyf'T'" ^ ^ 

^ ^ I ^minC xj €6 T ^minCG 1 €8P(Ti OO). 

\ V 4min / 

Noting that V(x(to)) < oc, it then follows from the Borel-Cantelli Lemma [23, Chapter III, 
Theorem 1.2] that P{3N < oo s.t. Vz > N, T'' < oo) = 0 [21, (2.28)-(2.31)]. This implies 
that P (limsup^^^ d{x{t),V) > 65 ) = 0 which contradicts (T39). Therefore ps = 0, that 
is 

P (lim inf d{x{t),'D) = 0 and lim sup d{x{t),T>{0)) > 0 ) =0. 

V t^oo J 

Finally, given that the events associated with pi, p 2 , and p 3 are nonintersecting, it follows 
that Pi = 1. Summarizing, we have shown that 


pi = P ( limsupd(x(t),P) =01=1 

P 2 — P fliminf d(x(t),P) > o) = 0 

V t^oo J 

P 3 = P (lim inf d{x{t), P) = 0 and limsupd(a:(t), P(0)) > 0 ) =0, 

V t^oo J 

and thus claim 2 of the theorem has been proven. We wish to reiterate the fact that this 
sketch barrows heavily from [21]. In most instances intermediate steps where stopping 
times are needed were glossed over. It might be worthwhile to revisit this analysis in the 
future, seeing as there does not seem to be much literature along this direction other than 
the references here in. □ 


For more detail on stopping times see [24, Part 2, Chapter II]. 
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Remark 3. The definition given for uniformly bounded in probability in this work is 
stronger than any other definitions that could be found in the literature and follows the 
classic definition, see Definition 2 and compare it to Definition 11. A trajectory is uni¬ 
formly bounded in probability, if for all r > 0 there exists a B(r) >0 such that ||x(to)|| < t 
implies 

P ^sup||a:(t)|| < B{r)J = 1 

for all t >to. In [21, Definition 2.2] the following definition is given, for all r > 0 and any 
e > 0 there exists a B(r, e) > 0 such that ||x(to)|| < r implies 

P [[sup||a:(t)|| < B{r)^ = 1 - e. 

In [50, §1.4] the definition of uniformly bounded is defined as follows, 

supP(||a:(t)II > ^ 0 as ^ oc. 

t 

Our definition of asymptotic attractivity follows that of Deng and Kristie [21] as it is 
allready as strong as the classic definition. Again however the often cited work by Khas- 
minskii [50] gives a much weaker definition. Compare the following definitions for asymp¬ 
totically attracting 


P flim ||a:(i)|| = o) = 1 

Vt^oo / 

([21, Definition 2.2]) 

1™ P ( lira ||x(i)|| = o) = 1. 

CCQ^O Vt—>-oo / 

([50, Equation (5.15)]) 


It is imperative when discussing stability that uniform bounds are achieved. 

4.5 Diagonal Stability of Random Matrices 

Theorem 17. If An G is chosen as 

[An]ij ~ ^=====A/'(0,1), i ^ j 

V (2 + 

for any (5 > 0, and [Anju = — 1, then An is asymptotically almost surely diagonally stable. 

Proof. Consider the random matrix An G defined as 

[An]zj ~ . =A/'(0,1), i A 

V(2 + (5)n 

and [An]ii = 0. Then it follows from (T2) that Var[An]zj = ( 2 +s)n ^ A 3- Let 

Bn — AnP Aff^ then it follows that Var[Bn]ij = { 2 + 5 )n ^Len i A 3 0 otherwise. From 

Theorem 3 we have that 


— I ‘In 

sup Xi{Bn) = 2J (1 + o(l)) 

l<i<n y (2 + ojn 

< 2 ( 1 + 0 ( 1 )) 

as n ^ oo asymptotically almost surely. Noting that Aj^ + An = Bn — ‘2Inxn it follows 
that An P An < 0 asymptotically almost surely. □ 

Corollary 18. If An G i is chosen as 


[Anhjr^AfiO^ 


1 


{2pS)n 




for any 6 > 0, and [An]zz 


1 , then An is asymptotically almost surely diagonally stable. 
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5 Clustering Analysis and Ordination Methods 

In the following section clustering techniques are explored in detail as well as ordination 
techniques for visualizing data. 


5.1 Distance and Metrics 

Let X = [xi, X2, .. •, Xn]^ G then the Euclidian norm is defined as 


lkll = 



For a row vector the euclidian norm is similary defined. If no subscript is given with ll'H 
then we assume the Euclidian norm is being used. The following is the Euclidian distance 
function d{x, y) = ||x — y||. 

A common distance metric used in ecology is the Jensen-Shannon Distanee (JSD) 
metric [46,54,59] 

JSD(y,2) 4 (T40) 

where the Jensen-Shanon divergence, JS(y,z) = KL(y, + z)) + KL(z, + y)), is 
simply the symetrized version of the Kulback-Liebler directed divergence KL(y, z) = 
^zlog—. So as to not divide by zero, a pseudo count of le-10 is added to zero 
elements before performing the JSD. 


5.2 /c-Medoids 


Assume that one has a collection of samples X G where n is the total number of 

samples and p is the dimension of each sample. Let Xi G z = 1, 2,..., n, be the row 

vectors of X as defined below 


A = 


"Ar 

A2 


A 


n. 


(T41) 


We would like each Xj to belong to a unique cluster within a collection of k possible 
clusters Ci,C 2 ,...,The unique cluster that contains sample Xj is denoted C{j; k). If 
sample Xj is contained in cluster 3 then C{j) = C 3 . If one is given an a-priori number of 
clusters then it is possible to perform this task using the popular paradigm of /c-medoids. 
The paradigm works as follows. Initially, k samples are chosen at random as representative 
medoids. k clusters are then constructed by associating other samples to the nearest 
medoid. Within each cluster all elements are tested so as to see if a different sample has 
a smaller within cluster sum of distances. The element with the smallest within cluster 
sum of distances is chosen as the new medoid for that cluster. This process is performed 
for each cluster. New clusters are constructed with the k new medoids and the algorithm 
repeats again. The MATLAB command kmedoids performs /c-medoids with a variety of 
different methods depending on the number of samples. When the number of samples is 
less than 3000 MATLAB implements /c-medoids using the Partitioning Around Medoids 
(PAM) algorithm [49]. 


5.3 Silhouette Value 

For each sample Xj there is a corresponding silhouette value Sj G [—1,1] which is defined 
as follows 

A - ai(fc) 

^ max{aj{k),bj{k)} 

where aj is the average dissimilarity between sample j and all other samples within its 
own cluster, hj is the average dissimilarity between Xj and the elements of the nearest 
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cluster, and k is the total number of apriori designated clusters. These two quantities are 
now formally defined as follows 

' Xiecu-,k) 


Note that we have used \C{j) \ to denote the cardinality of C(j), the total number of samples 
in C{j). We similarly define the bj as 


bj{k) 


A 


mm —— 

Cm^C{j-,k) \Crr 


XiECm 


For a given sample set and a given number of clusters k > 2 there is a corresponding Sj{k). 
The Silhouette Index for a sample data set is then the maximum of the mean silhouette 
value for each total number of clusters. 

1 "" 

SI(X) ^max-y]sj(fc) 

i=i 


The optimal number of clusters is then argmax^ ^ Silhouette analysis is 

performed in MATALAB using the command silhouette. 


5.4 Variance Ratio Criterion and the Calihski-Harabasz Index 


We now define the Variance Ratio Criterion (VRC) which holds for any distance function. 
When the Euclidian metric is used the VRC is referred to as the Calinski-Harahasz (CH) 
index [12]. As before, assume that a collection of samples has already been grouped into 
k clusters. The VRC is defined as 


YRC{k) = 


BG{k) n-k 
WG{k) k - 1 


where BG is the Between Group variance and WG is the Within Group variance defined 
below. 


BG(fc) ^ ^ 


1 l<-il |Bm| XieCj Xiec. 

m>j 


Y Y d{Xi,Xef 


wG(fc) ^ 


n n 


-3TT E E d{X,,Xe)\ 


^|Cil(|Cj| 1) ^ 

X^^Cj 


5.5 Principle Coordinate Analysis 


The purpose of Principle Goordinates Analysis (PCoA) is to represent a collection of high 
dimensional data in a lower dimension. Once again assume that one has a collection of 
samples X G where n is total number of samples and p is the dimension of each 

sample. Let Xi G R^^^, z = 1, 2,... ,n be the row vectors of X as defined in (T41). The 
question answered in this section is how one obtains ayGR’^^^,A:<n with Yi G R^^^, 
z = 1, 2 ,..., n defined as follows 

r^ii 


Y = 


^2 


L>"nJ 

that is a faithful representation of X. We begin by defining the dissimilarity between sam¬ 
ples z and j as Sij = d{Xi, Xj). Then the goal of this method is to find Y such that d(V, Yj) 
is similar to d{Xi,Xj) for the distance measure of interest. Let D = and 


B={lnxn~n ^l„ln)D{I„x„ ~ n 








6. MODELING 


61 


where In is an n-dimensional column vector with each entry equal to 1. The n x k 
dimensional representation of the n x p sample data is then [82, Chapter 5] 

Y = [qi a/a7, q 2 ■ ■ ■, qk \/^] • 

where B = QAQ~^ is the eigenvalue decomposition with eigenvalues Xi G M, and normal¬ 
ized eigenvectors qi G 'ML for i = 1, 2,..., n with Q = [^i, ^ 2 , • • •, qn] and the 

diagonal eigenvalue matrix. Due to the fact that B is symmetric all eigenvalues and eigen¬ 
vectors will be real valued. It is furthermore assumed that the eigenvalues are arranged 
such that Ai > A 2 > • • • > An- When the Euclidian distance function is used B is always 
positive semi-definite. However, it may be possible to have negative eigenvalues when 
other distance functions are used and thus it is necessary to check and make sure that 
that Xi for each i G {1, 2,..., /c} is positive for the /c-dimensional space of interest. The 
command that performs this task in MATLAB is cmdscale. MATLAB also rotates the 
data, but that has no affect on the dissimilarity measure, as rotation is distance preserving 
in Euclidian space. 

5.6 Principle Component Analysis 

Principle Component Analysis (PCA) is similar in spirit to PCoA in that one wishes to 
represent a set of data in a lower dimension. Assume that one has a collection of samples 
X G where n is total number of samples and p is the dimension of each sample. As 

before the goal is to find a T G k < p that is representative of the original data. 

Assume without loss of generality that the columns of X have mean zero. Let 

be the singular value decomposition of X where E G R’^^^ is a diagonal matrix with the 
singular values cr^ = [EJ^* where ai > aj for i < j. The columns of G are the left 

singular vectors and the columns of IT G R^^^ are the right singular vectors. The reduced 
order representation is then defined as 

T = XWl:n,l:k 

where Wi-.np-.k contains only the first k columns of IT. Note that Wi:n,i:k can be thought 
of as a right hand side projection operator and can be used to project subsequent mea¬ 
surements into pre-existing principle components. Also, note that when the Euclidean 
distance is used in PCoA it is equivalent to PC A. The command that performs this task 
in MATLAB is pea. 

6 Modeling 

Two approaches for modeling individual human microbial communities are now presented. 
The first model is the same as that presented in the main text and from this point forward 
is referred to as the universal model. This model assumes that all species interact in a 
universal manner independent of the host. The second model allows for there to be multiple 
possible interaction strengths and growth rates for species. This model will be referred 
to as the multiple set model We use Eigure T4 as a visual reference when discussing the 
modeling paradigms. 

6.1 Universal Model 

Consider a universal species pool indexed by a set of integers S = {1, ..., n}, a global 
n X n matrix A representing all possible pairwise interactions between species, and a 
universal vector r of size n containing the growth rates for all the n species. The global 
variables for our ecological system are completely defined by the triple (S,A,r). Then 
q Local Communities (LCs) are defined by sets which are subsets of S denoting the 
specific microbes present in LCj^, u = 1,...,^. Eor simplicity we assume that each LC 
contains only p species (p < n), randomly selected from the universal species pool. The 
GLV dynamics for each LC is given by 

LCn : = diag x^'^^t)^ 


(T42) 
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where the LC specific interaction matrix and growth vector are defined as A^^'^ = 
and , respectively. That is, A^^^ is obtained from A by only taking the rows 

and columns of A that are contained in the set A similar procedure is performed 

in order to obtain Finally there is a map rrij^ that takes the abundances in the 
index and carries them to the universal index in S giving the vector = rrijy 
which results in for j = 1 ,... ,p and = 0 for j in S \ . This modeling 

procedure is inspired by [18]. A toy example of how this model is used to construct a 
single local community is now given 

Example 1. Consider a global set of 4 species and thus S = {1, 2, 3,4}. Let 
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Universal Dynamics 



Multiple Sets of Dynamics 



LCi LC2 LCs LCg 

Supplementary Text Figure T4: Modeling Paradigms. Two different modeling 
paradigms are presented. In the universal paradigm there is a universal list of 
species S, a universal interaction matrix A, containing all the pairwise interaction 
strengths, and a universal growth rate vector r. Then the dynamics of each Local 
Community LC is determined by the collection of species in that community. In 
the multiple set paradigm the dynamics are not only determined by the collection 
of species present but the dynamics are not universal and come from a collection of 
possible dynamics. In the example above LCi dynamics are derived from the first 
dynamic pair (Ai,ri), LC 2 dynamics are derived from (A 2 ,r 2 ), LC 3 dynamics 
are derived from (A^jF^), and LCg dynamics are derived from the dynamic pair 
(Ai, ri). 
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Then, if the first local community has the species list S^^'^ = {1,3,4} then it follows that 


■ -1 

0.4 

-o.r 

and = 
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6.2 Mutiple Set Model 

As before there is a universal set of species S = {l,...,n}, but now there are i pairs 
of possible global dynamics {(Ai, ri), (A 2 , r 2 ),..., (A^, r^)}. The q LCs are defined by 
sets S^^'^ , u = 1 ,... ,q and a map w : {1, 2 ,..., q} ^ {1, 2 ,..., ^} which determines which 
model pair the LC is derived from. The GLV dynamics for each LC are given by (T42) 
where and , respectively. That is, is obtained 

from by only taking the rows and columns of that are contained in the set 

A similar procedure is performed in order to obtain The map rriu that takes the 
abundances in the index and carries them to the universal index S is defined as 
before = mu{x^^^). It is worth noting that results from this modeling paradigm are 
not presented in the main text, because the emergence of community types in this case is 
trivial (see §7.1) [52]. 

7 Simulation Results and Analysis 

We are finally in a position to address the main topic of this work. Revealing properties 
that allow for community types to arise. Four case studies are now presented. The first case 
study will use the multiple set model paradigm in Seciton 6 . 2 , and it will be shown that 
under this paradigm that clustering of the steady states occurs trivially. The subsequent 
three case studies use the universal model from Section 6.1. The case studies explore 
heterogeneity (in terms of interaction strength and network topology), mean degree of the 
network, and community overlap and how these affect the existence of community types 
for our dynamics of interest. 

7.1 Multiple Set Models 

For the multiple set model study three universal pairs {(Ai, ri), (A 2 , r 2 ), (A 3 , ra)} were 
considered. For this study there was a total of 100 global species. Each Ai was a 100 by 
100 matrix with values drawn from a normal distribution with mean 0 and variance 0.0049. 
Then the diagonal elements of each matrix are set to -1, that is [Ai]jj = —1, i = 1,2,3, 
j = 1,2, ...,100. Note that with the above selection, the variance of the off diagonal 
elements of Ai satisfy the following bound 0.0049 < ^Too’ which from Corollary 18 is 
the sufficient condition on the variance for asymptotic almost sure diagonal stability when 
n — 100. Also, note the interaction network is implicitly a complete graph without any 
structural heterogeneity. Finally, each had elements drawn from a uniform distribution 
between 0 and 1 . 

From the three pairs above 600 local communities were constructed, 200 from each 
universal pair. Each of the local communities contained 80 species. The dynamics were 
then simulated for 100 seconds with initial conditions drawn from Z//(0,1). The abundances 
of the species in the communities were then normalized, and the relative abundances of 
the 600 local communities were clustered using /c-medoids from Section 5.2 and silhouette 
indexed as defined in Section 5.3, each using the Jensen-Shannon metric in (T40). A 
principle coordinate plot of the results is given in Eigure T5 with the elements colored 
according to cluster assignment from /c-medoids. It is clear that there are three clusters. 
Eurthermore each cluster exactly coincides with the universal A^ that the local community 
dynamics were obtained from. 

7.2 Universal Model 

Eor the three universal case studies the global interaction matrix is defined as 

A = NH o Gs, 


which contains four components. 
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(i) N G is the nominal component where each element is sampled from either a 

normal distribution or a uniform distribution. 

(ii) The matrix H is a diagonal matrix that captures the overall interaction heterogeneity 
of different species. When interaction strength heterogeneity is employed the diago¬ 
nal elements of H are drawn from a power-law distribution (T3) with exponent —a, 
[H.]ii ~ P{oi), which are subsequently normalized so that the mean of the diagonal 
is equal to 1. Without interaction heterogeneity H is simply the identity matrix. 

(iii) The matrix G is the adjacency matrix of the underlying ecological network: [G]ij = 1 
if species i is affected by the presence of species j and 0 otherwise. When the 
underlying network is a complete digraph all elements in G are equal to 1. For 
details on the construction of G when the underlying network is Erdds-Renyi or a 
power-law digraph see Sections 3.2.1 and 3.2.2 respectively. 

(iv) The last element s is simply a scaling factor between 0 and 1. 

As before we set = — 1 to ensure one of the necessary conditions for diagonal stability 
of a matrix is satisfied (Theorem 10 and Corollary 11). Finally, the elements in the global 
growth rate vector are dehned as 

r = h o n 

where n is the nominal component taken from a uniform distribution, and h captures the 
growth rate heterogeneity. When there is growth rate heterogeneity h is drawn from a 
power-law distribution with exponent —a and subsequently normalized to have a mean of 
1. Without growth rate heterogeneity h is simply a column vector of ones. Note that h is 
not included in the main text, and is only used in one scenario of one case study within 
this section. 

7.2.1 Universal Model: Heterogeneity Study 

Next is a systematic study of heterogeneity and its effects on the dynamics in the universal 
model. Eight different scenarios were tested in this study. For each of the following 
scenarios the number of global species is n = 100. Table T1 outlines the differences 
between the scenarios. 


clusters=3, Si=0.39 



Principal Coordinate 1 


Supplementary Text Figure T5: Principle coordinates for multiple model case 
study. 
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Scenario 1. [N]ij ~ 0.5,0.5), [H]ii ~ where a G [1.2,7] and subsequently 

normalized to have a mean of 1, G is the adjacency matrix for a complete 
digraph and thus all entries are equal to 1, and the scaling factor is set 
to s = 0.07. There is no growth rate heterogeneity and thus the column 
vector h is all ones and finally [n]^ ~ Z//(0,1). 

Scenario 2. The same as Scenario 1 but with [Ni]ij ~ A/’(0,1) and s = 0.07. 

Scenarios. [Ni]ij ~ 0.5,0.5), [H]^^ ~ ”P(<^) where a G [1.2,7] and subsequently 

normalized to have a mean of 1, G is the adjacency matrix for an Erdds- 
Renyi digraph with a mean out-degree of 10, and the scaling factor is set to 
s = 0.5. There is no growth rate heterogeneity and thus the column vector 
h is all ones and finally [n]^ ~ ^7(0,1). For more details on Erdds-Renyi 
digraphs see §3.2.1. 


Scenario 4. The same as Scenario 3 but with [Ni]ij ~ A/’(0,1) and s = 0.1. 

Scenario 5. [N]ij ~ A/'(0,1), H is the identity matrix, G is the adjacency matrix for 
a digraph with the out-degree drawn from a power-law distribution V{a) 
with a mean out-degree of 10 where a G [1.2,7], and the scaling factor is 
set to s = 0.2. There is no growth rate heterogeneity and thus the column 
vector h is all ones and [n]i ~ Z//(0,1). For more details on power-law 
out-degree digraphs see §3.2.2. 


Scenario 6. [N]^^ ~ A/’(0,1), [H]^^ ~ where a G [1.2, 7] and subsequently normal¬ 
ized to have a mean of 1, G is the adjacency matrix for a digraph with 
the out-degree drawn from a power-law distribution V{a) with a mean 
out-degree of 10 where a G [1.2, 7], and the scaling factor is set to s = 0.1. 
The power-law degree distribution and interaction strength heterogeneity 
are uncoupled. There is no growth rate heterogeneity and thus the column 
vector h is all ones and [n]i ~ 14(00)- Foi* more details on power-law 
out-degree digraphs see §3.2.2. 


Scenario 7. [N]^^ ~ A/’(0,1), [H]^^ ~ ”?^(<^) where a G [1.2,7] and subsequently nor¬ 
malized to have a mean of 1, G is the adjacency matrix for a digraph 
with the out-degree drawn from a power-law distribution V(a) with a 
mean out-degree of 10 where a G [1.2,7], and the scaling factor is set 
to s = 0.02((a + 1). The power-law degree distribution and interaction 
strength heterogeneity are coupled so that the node with the highest out- 
degree is also the node with the largest interaction strength scaling. There 
is no growth rate heterogeneity and thus the column vector h is all ones 
and [n]i ~ 14(0, 1). For more details on power-law out-degree digraphs see 
§3.2.2. 


Scenario 8. [N]^^ ~ A/'(0,1), H is the identity matrix, G is the adjacency matrix for 
an Erdds-Renyi digraph with a mean out-degree of 10, and the scaling 
factor is set to s = 0.1. There is growth rate heterogeneity and thus the 
column vector h has elements drawn from a power-law distribution V(a) 
and subsequently normalized to have a mean of 1 and [n]^ ~ 14(0, 1). For 
more details on Erdds-Renyi digraphs see §3.2.1. 


For each of the eight scenarios above and for every a, q — 500 local communities were 
generated each with p = 80 species selected at random from the n = 100 global species. 
The dynamics as described above following the modeling paradigm in Section 6.1 were 
then simulated for 100 seconds with initial conditions drawn from U(0, 1). If any of the 
500 simulations crashed due to instability or if the norm of the terminal discrete time 
derivative was greater than 0.01 then that local community was excluded from the rest 
of the study. Those simulations that finished without crashing and with small terminal 
discrete time derivative were deemed steady. Less than 1% of simulations were deemed 
unstable. The abundances of the species in the communities were then normalized, and 
the relative abundances of the 500 local communities were clustered using /c-medoids from 
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Section 5.2 and silhouette indexed as defined in Section 5.3, each using the Jensen-Shannon 
distance metric in (T40). 

The results from the above scenarios are presented in Figures T6-T13. The first plot is a 
comprehensive clustering analysis of the steady state values obtained from the simulations. 
The x-axis denotes the heterogeneity value a. The box plots are the silhouette values 
pertaining to the number of clusters for which the silhouette index was defined. The total 
number of clusters pertaining to the silhouette index is denoted on the top x-axis. The 
second row is a principle coordinate analysis of the steady state values obtained at three 
different heterogeneity values. Clusters are color coded to match the optimal clustering 
from /c-medoids. The third row of the figure plots 


max real 




as a function of a. The fourth row is a plot of Xi{A^^^) for i G {1,2,..., 80} and j G 
{1,2,..., 500} at three values of a. 

Figures T6 and T7 show that regardless of whether the nominal component is drawn 
from a uniform or a normal distribution the increase in interaction strength heterogeneity 
(decreasing a) leads to steady state clustering in the data. Figures T8 and T9 illustrate 
that the same phenomenon also holds for Erdos-Renyi digraphs as well. When the under¬ 
lying graph topology follows a power-law degree distribution, but there is no interaction 
strength heterogeneity, clustering of steady states is not observed, see Figure TIO. Figure 
Til illustrates the fact that when there is interaction strength heterogeneity and network 
degree heterogeneity it is possible to have clustering of steady states when a is in the range 
[3,1]. However the trend is not smooth and is inconsistent, as compared to Figures T6-T9. 
This is due to the fact that when the underlying interaction network is being constructed, 
there is no guarantee that the high-degree node will also have a large interaction strength. 
For instance, if a node with no out edges is randomly selected to have high interaction 
strength scaling, then the impact of that node on the rest of the nodes is still zero. When 
the interaction strength heterogeneity is coupled with the out-degree for a power-law out- 
degree digraph then the trends from Figures T6-T9 are recovered. Finally, if the growth 
rates of the species are derived from a power-law distribution, then clustering of steady 
states also occurs as a decreases. 

Rows three and four illustrate the spectrum of A^^\ and are included so that we can 
infer the asymptotic stability of the system for certain paradigms. Note that regardless 
of whether the nominal interactions are drawn from a normal distribution or a uniform 
distribution when a is large the spectrum represents a uniform disk in the complex plain as 
predicted by Theorem 1. Furthermore, for scenario 2 with s = 0.07 and [N]ij ~ A/’(0,1) it 
follows that Var[A]ij < for n = 100. From Theorem 17 it follows that A is diagonally 
stable, in a probabilistic sense. Then invoking Theorem 10 we know that any principle 
minor of A is diagonally stable as well. Therefore, in Scenario 2 for large a each A^^'^ is 
diagonally stable, in a probabilistic sense. 

Also, for all of the scenarios with interaction heterogeneity all of the eigenvalues of A, 
and consequently A^^\ converge to —1 as a tends to 1. In the limit of a tending to 1, only 
one of the columns of A has non-zero values off the diagonal. Therefore, in the limit of a 
tending to 1 the following inequality holds A^P + PA < 0 where P = [1, e,..., e]^ and e is 
sufficiently small. ^ Thus A is diagonally stable in the limit of a tending to 1. Therefore, 
for low interaction strength heterogeneity and for high interaction strength heterogeneity 
one of the necessary conditions for uniform asymptotic stability in the positive orthant is 
satisfied, see Theorem 5. Note that even when the A^^'^ are not diagonally stable it does 
not imply that the state trajectory is unstable. 


7.2.2 Universal Model: Sparsity Study 

Next is a study where the mean degree of the Erdos-Renyi digraph along with the in¬ 
teraction strength heterogeneity is varied. Table T2 outlines the differences between the 
scenarios. The details of the study are as follows: [N]ij ~ 0.5, 0.5), [H]ii ~ ^{a) 

where a G [1.2,7] and subsequently normalized to have a mean of 1, G is the adjacency 

^Without loss of generality we have assumed that the first column of A is the highly weighted column. 
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matrix for an Erdos-Renyi digraph with a mean out-degree of 
d e {1, 3, 5, 7, 9, 11, 13, 15, 17, 19}, 

and the scaling factor is set to s = I/a/S. There is no growth heterogeneity and thus the 
column vector h is all ones and finally [n]^ ~ IL{0, 1). For each of the ten scenarios above 
and for every a the same procedure as in §7.2.1 was carried out and the results are shown 
in Figures T14-T23. 

From Figures T14-T23 it can be concluded that so long as the mean out-degree of 
the ER digraph is greater than 2 the steady states of the GLV model increase in SI as a 
decreases. That is, the same trends as observed in the previous study hold, so long as the 
ER digraph is connected. When the mean degree is 1 for an ER graph it is very unlikely 
that the graph will be connected [29]. Thus for a digraph with mean out-degree 1, it is 
even more unlikely that it will be connected. When the underlying digraph G has many 
isolated nodes then the scaling of interaction strengths does not influence as many other 
species in the GLV system and thus clustering is not observed. 

7.2.3 Universal Model: Community Size Study 

In all previous studies each local community contained 80 species. For this study the size 
of the local communities take on values in the following set 

p e {100, 99, 95, 90, 80, 70, 60, 50}. 

Details of the study are as follows: [N]ij ~ 0.5, 0.5), [H]^^ ~ where a G [1.2,7] 

and subsequently normalized to have a mean of 1, G is the adjacency matrix for an Erdos- 
Renyi digraph with a mean out-degree of 10, and the scaling factor is set to s = 0.5. 
There is no growth heterogeneity and thus the column vector h is all ones and finally 
[n]i ~ ^7(0,1). Table T3 outlines the differences between the scenarios. For each of the 
eight scenarios above and for every a the same procedure as in §7.2.1 was carried out and 
the results are shown in Figures T24-T31. 

These results show that the trends observed in the earlier studies still hold when the 
community sizes are varied between 95 and 50 species. Figures T26-T31. As the community 
sizes approach 50 the trend of increased clustering with increased heterogeneity is less 
significant. When the number of species in the LGs approaches the number of species 
in the met a-community, 100, the results do not follow the same trends as before and the 
Silhouette Indices are near 1, independent of the interaction strength heterogeneity. In 
Figure T24 all the simulations are identical, each LG has the same 100 species but with 
different initial conditions. Due to asymptotic stability all of the simulations converge to 
the same steady state (only those simulations for a between 2.2 and 3 have the potential 
to be unstable). Even though the Silhouette Indices are near 1 across the heterogeneity 
spectrum, all of the steady state values are within 10~^ to 10~^ in terms of the first 
two principle coordinates. When all of the LGs differ by only one species. Figure T25, 
once again the Silhouette Indices are near 1, yet the PGoA illustrates that most of the 
samples are very close together with just a few outliers. Figures T24 and T25 illustrate the 
sometimes confounding results when clustering analysis and PGoA are performed [44,67]. 
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Supplementary Text Figure T6: Universal Model Heterogeneity Study Scenario 
1 in Table Tl. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxj^j real(Ai(A['^])) as a function of a. The fourth row is a plot of Xi{A^^^)for 
i G {1, 2, . . . , 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T7: Universal Model Heterogeneity Study Scenario 
2 in Table Tl. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxj^j real(Ai(A['^])) as a function of a. The fourth row is a plot of Xi{A^^^)for 
i G {1, 2, . . . , 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T8: Universal Model Heterogeneity Study Scenario 
3 in Table Tl. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxj^j real(Ai(A['^])) as a function of a. The fourth row is a plot of Xi{A^^^)for 
i G {1, 2, . . . , 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T9: Universal Model Heterogeneity Study Scenario 
4 in Table Tl. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a:-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top cc-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /c-medoids. The third row of the figure 
plots maxij real(Ai(A['^])) as a function of a. The fourth row is a plot of Ai(At-^^)for 
i G {1, 2, . . . , 80} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure TIO: Universal Model Heterogeneity Study Sce¬ 
nario 5 in Table Tl. The first plot is a comprehensive clustering analysis of the 
steady state values obtained from the Lotka-Volterra simulations. The a^-axis de¬ 
notes the heterogeneity value a. The box plots are the silhouette values pertaining 
to the number of clusters for which the silhouette index (maximum over the num¬ 
ber of clusters of the mean silhouette value for each given total number of cluster) 
was defined. The total number of clusters pertaining to the silhouette index is de¬ 
noted on the top a:-axis. The second row is a principle coordinate analysis of the 
steady state values obtained at three different heterogeneity values. Clusters are 
color coded to match the optimal clustering from /c-medoids. The third row of the 
figure plots may^ij real(Ai)) as a function of a. The fourth row is a plot of 
Ai(A['^])for i G {1, 2, . . . , 80} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure Til: Universal Model Heterogeneity Study Sce¬ 
nario 6 in Table Tl. The first plot is a comprehensive clustering analysis of the 
steady state values obtained from the Lotka-Volterra simulations. The a^-axis de¬ 
notes the heterogeneity value a. The box plots are the silhouette values pertaining 
to the number of clusters for which the silhouette index (maximum over the num¬ 
ber of clusters of the mean silhouette value for each given total number of cluster) 
was defined. The total number of clusters pertaining to the silhouette index is de¬ 
noted on the top a:-axis. The second row is a principle coordinate analysis of the 
steady state values obtained at three different heterogeneity values. Clusters are 
color coded to match the optimal clustering from /c-medoids. The third row of the 
figure plots may^ij real(Ai)) as a function of a. The fourth row is a plot of 
Ai(A['^])for i G {1, 2, . . . , 80} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure T12: Universal Model Heterogeneity Study Sce¬ 
nario 7 in Table Tl. The first plot is a comprehensive clustering analysis of the 
steady state values obtained from the Lotka-Volterra simulations. The a:-axis de¬ 
notes the heterogeneity value a. The box plots are the silhouette values pertaining 
to the number of clusters for which the silhouette index (maximum over the num¬ 
ber of clusters of the mean silhouette value for each given total number of cluster) 
was defined. The total number of clusters pertaining to the silhouette index is de¬ 
noted on the top cc-axis. The second row is a principle coordinate analysis of the 
steady state values obtained at three different heterogeneity values. Clusters are 
color coded to match the optimal clustering from /c-medoids. The third row of the 
figure plots mayiij real(Ai)) as a function of a. The fourth row is a plot of 
Ai(A['^])for i G {1, 2, . . . , 80} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure T13: Universal Model Heterogeneity Study Sce¬ 
nario 8 in Table Tl. The first plot is a comprehensive clustering analysis of the 
steady state values obtained from the Lotka-Volterra simulations. The a^-axis de¬ 
notes the heterogeneity value a. The box plots are the silhouette values pertaining 
to the number of clusters for which the silhouette index (maximum over the num¬ 
ber of clusters of the mean silhouette value for each given total number of cluster) 
was defined. The total number of clusters pertaining to the silhouette index is de¬ 
noted on the top a:-axis. The second row is a principle coordinate analysis of the 
steady state values obtained at three different heterogeneity values. Clusters are 
color coded to match the optimal clustering from /c-medoids. The third row of the 
figure plots maxij real(Ai)) as a function of a. The fourth row is a plot of 
Ai(A['^])for i G {1, 2, . . . , 80} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure T14: Universal Model Sparsity Study Scenario 1 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A['^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T15: Universal Model Sparsity Study Scenario 2 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A[-^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T16: Universal Model Sparsity Study Scenario 3 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A[-^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T17: Universal Model Sparsity Study Scenario 4 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A['^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T18: Universal Model Sparsity Study Scenario 5 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A['^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T19: Universal Model Sparsity Study Scenario 6 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A[-^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T20: Universal Model Sparsity Study Scenario 7 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A[-^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T21: Universal Model Sparsity Study Scenario 8 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A[-^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T22: Universal Model Sparsity Study Scenario 9 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A[-^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T23: Universal Model Sparsity Study Scenario 10 
in Table T2. The first plot is a comprehensive clustering analysis of the steady 
state values obtained from the Lotka-Volterra simulations. The a^-axis denotes the 
heterogeneity value a. The box plots are the silhouette values pertaining to the 
number of clusters for which the silhouette index (maximum over the number of 
clusters of the mean silhouette value for each given total number of cluster) was 
defined. The total number of clusters pertaining to the silhouette index is denoted 
on the top a:-axis. The second row is a principle coordinate analysis of the steady 
state values obtained at three different heterogeneity values. Clusters are color 
coded to match the optimal clustering from /e-medoids. The third row of the figure 
plots maxij real(Ai(A[-^^)) as a function of a. The fourth row is a plot of 
for i G {1,2,..., 80} and j G {1,2,..., 500} at three values of a. 
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Supplementary Text Figure T24: Universal Model Community Size Overlap 
Study Scenario 1 in Table T3. The first plot is a comprehensive clustering analysis 
of the steady state values obtained from the Lotka-Volterra simulations. The x- 
axis denotes the heterogeneity value a. The box plots are the silhouette values 
pertaining to the number of clusters for which the silhouette index (maximum over 
the number of clusters of the mean silhouette value for each given total number of 
cluster) was defined. The total number of clusters pertaining to the silhouette index 
is denoted on the top a:-axis. The second row is a principle coordinate analysis of 
the steady state values obtained at three different heterogeneity values. Clusters 
are color coded to match the optimal clustering from /c-medoids. The third row of 
the figure plots may^ij real(Ai(Al[-^])) as a function of a. The fourth row is a plot of 
) for i G {1, 2, . . . , p} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure T25: Universal Model Community Size Overlap 
Study Scenario 2 in Table T3. The first plot is a comprehensive clustering analysis 
of the steady state values obtained from the Lotka-Volterra simulations. The x- 
axis denotes the heterogeneity value a. The box plots are the silhouette values 
pertaining to the number of clusters for which the silhouette index (maximum over 
the number of clusters of the mean silhouette value for each given total number of 
cluster) was defined. The total number of clusters pertaining to the silhouette index 
is denoted on the top a:-axis. The second row is a principle coordinate analysis of 
the steady state values obtained at three different heterogeneity values. Clusters 
are color coded to match the optimal clustering from /c-medoids. The third row of 
the figure plots max^^j real(Ai(A['^^)) as a function of a. The fourth row is a plot of 
for i G {1, 2, . . . ,p} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure T26: Universal Model Community Size Overlap 
Study Scenario 3 in Table T3. The first plot is a comprehensive clustering analysis 
of the steady state values obtained from the Lotka-Volterra simulations. The x- 
axis denotes the heterogeneity value a. The box plots are the silhouette values 
pertaining to the number of clusters for which the silhouette index (maximum over 
the number of clusters of the mean silhouette value for each given total number of 
cluster) was defined. The total number of clusters pertaining to the silhouette index 
is denoted on the top a:-axis. The second row is a principle coordinate analysis of 
the steady state values obtained at three different heterogeneity values. Clusters 
are color coded to match the optimal clustering from /c-medoids. The third row of 
the figure plots max^^j real(Ai(Al['^^)) as a function of a. The fourth row is a plot of 
for i G {1, 2, . . . ,p} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure T27: Universal Model Community Size Overlap 
Study Scenario 4 in Table T3. The first plot is a comprehensive clustering analysis 
of the steady state values obtained from the Lotka-Volterra simulations. The x- 
axis denotes the heterogeneity value a. The box plots are the silhouette values 
pertaining to the number of clusters for which the silhouette index (maximum over 
the number of clusters of the mean silhouette value for each given total number of 
cluster) was defined. The total number of clusters pertaining to the silhouette index 
is denoted on the top a:-axis. The second row is a principle coordinate analysis of 
the steady state values obtained at three different heterogeneity values. Clusters 
are color coded to match the optimal clustering from /c-medoids. The third row of 
the figure plots max^^j real(Ai(A['^^)) as a function of a. The fourth row is a plot of 
for i G {1, 2, . . . ,p} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure T28: Universal Model Community Size Overlap 
Study Scenario 5 in Table T3. The first plot is a comprehensive clustering analysis 
of the steady state values obtained from the Lotka-Volterra simulations. The x- 
axis denotes the heterogeneity value a. The box plots are the silhouette values 
pertaining to the number of clusters for which the silhouette index (maximum over 
the number of clusters of the mean silhouette value for each given total number of 
cluster) was defined. The total number of clusters pertaining to the silhouette index 
is denoted on the top a:-axis. The second row is a principle coordinate analysis of 
the steady state values obtained at three different heterogeneity values. Clusters 
are color coded to match the optimal clustering from /c-medoids. The third row of 
the figure plots max^^j real(Ai(Al['^^)) as a function of a. The fourth row is a plot of 
for i G {1, 2, . . . ,p} and j G {1, 2, . . . , 500} at three values of a. 












Imaginary max* real(Ai) Principal Coordinate 2 Silhouette Values 


94 


CHAPTER 5. SUPPLEMENTARY TEXT 


Number of Clusters 

222222322223322232222222222222 




q:= 1.2, clusters=2, SI=0.55 

0.6 -.-^- 


CM 



Principal Coordinate 1 


Principal Coordinate 1 


Principal Coordinate 1 


“I-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1—r 


0.5 - 


* * 


*^******* 


* * * ¥ 


* * * 


-0.5 - 


* * ^ 


I I I I I I I I I I I I I I I I I I I I I I I I I I I I ^ 

7 6.8 6.6 6.4 6.2 6 5.8 5.6 5.4 5.2 5 4.8 4.6 4.4 4.2 4 3.8 3.6 3.4 3.2 3 2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 


1 

0.5 

0 

■0.5 

-1 

-2 -1 0-2 -1 0 
Real Real 


a=7 


a=2.6 



a=1.2 



Supplementary Text Figure T29: Universal Model Community Size Overlap 
Study Scenario 6 in Table T3. The first plot is a comprehensive clustering analysis 
of the steady state values obtained from the Lotka-Volterra simulations. The x- 
axis denotes the heterogeneity value a. The box plots are the silhouette values 
pertaining to the number of clusters for which the silhouette index (maximum over 
the number of clusters of the mean silhouette value for each given total number of 
cluster) was defined. The total number of clusters pertaining to the silhouette index 
is denoted on the top a:-axis. The second row is a principle coordinate analysis of 
the steady state values obtained at three different heterogeneity values. Clusters 
are color coded to match the optimal clustering from /c-medoids. The third row of 
the figure plots max^^j real(Ai(Al['^^)) as a function of a. The fourth row is a plot of 
for i G {1, 2, . . . ,p} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure T30: Universal Model Community Size Overlap 
Study Scenario 7 in Table T3. The first plot is a comprehensive clustering analysis 
of the steady state values obtained from the Lotka-Volterra simulations. The x- 
axis denotes the heterogeneity value a. The box plots are the silhouette values 
pertaining to the number of clusters for which the silhouette index (maximum over 
the number of clusters of the mean silhouette value for each given total number of 
cluster) was defined. The total number of clusters pertaining to the silhouette index 
is denoted on the top a:-axis. The second row is a principle coordinate analysis of 
the steady state values obtained at three different heterogeneity values. Clusters 
are color coded to match the optimal clustering from /c-medoids. The third row of 
the figure plots max^^j real(Ai(Al['^^)) as a function of a. The fourth row is a plot of 
for i G {1, 2, . . . ,p} and j G {1, 2, . . . , 500} at three values of a. 
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Supplementary Text Figure T31: Universal Model Community Size Overlap 
Study Scenario 8 in Table T3. The first plot is a comprehensive clustering analysis 
of the steady state values obtained from the Lotka-Volterra simulations. The x- 
axis denotes the heterogeneity value a. The box plots are the silhouette values 
pertaining to the number of clusters for which the silhouette index (maximum over 
the number of clusters of the mean silhouette value for each given total number of 
cluster) was defined. The total number of clusters pertaining to the silhouette index 
is denoted on the top a:-axis. The second row is a principle coordinate analysis of 
the steady state values obtained at three different heterogeneity values. Clusters 
are color coded to match the optimal clustering from /c-medoids. The third row of 
the figure plots max^^j real(Ai(Al['^^)) as a function of a. The fourth row is a plot of 
for i G {1, 2, . . . ,p} and j G {1, 2, . . . , 500} at three values of a. 
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